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Semiconductor  Measurement  Technology : 
DISTRIB  I,  An  Impurity  Redistribution  Computer  Progrcun 

by 

David  Gilsinn  and  Richard  Kraft 

This  report  provides  docxjmentation  of  a  computer  program 
which  calculates  the  redistribution  of  impurities  in 
silicon  during  a  single  oxidation  step.     The  dociomenta- 
tion  provides:     (1)  a  physical  and  mathematical  descrip- 
tion of  the  redistribution  process,    (2)  a  detailed  de- 
scription of  the  discretization  of  the  appropriate  par- 
tial differential  equations,  and  (3)  a  complete  descrip- 
tion of  the  FORTRAN  program  for  computing  the  solution. 

Key  Words:     Diffusion;   electronic  technology;  impurity 
distribution;  material  transport;   segregation;  semi- 
conductor technology. 

1 .  INTRODUCTION 

1.1.  Main  Objectives  of  the  Documentation 

The  principal  objectives  of  this  documentation  are: 

(1)  to  give  a  detailed  derivation  of  the  algebraic  equations  that 
are  used  in  the  computer  program, 

(2)  to  give  in  detail,  a  step-by-step  description  of  the  way  the 
computation  in  each  computational  block  of  the  progreun  is  car- 
ried^ out,  and 

(3)  to  provide  a  simple  example  illustrating  how  to  use  the  pro- 
gram. 

The  format  of  this  documentation  is  modeled  after  the  partial  differential 
equation  computer  prograjn  format  description  in  [9]  in  the  FORTRAN  IV  lan- 
guage . 

In  fabrication  processes  redistribution  of  impurities  takes  place  in  mul- 
tiple oxidation  and  diffusion  steps   (each  step  being  characterized  by  con- 
stant process  parameters) .     DISTRIB  is  aimed  at  simulating  this  complete 
process;  however,  due  to  the  complexity  of  the  total  problem,  it  was  de- 
cided to  only  document  a  fundamental  oxidation  step  here.     This  ejcplains 
the  appellation  DISTRIB  VERSION  1. 

The  next  section  introduces  the  user  to  DISTRIB  VERSION  1. 

1.2.  Tips  on  Where  to  Locate  the  Key  Parts  of  This  Documentation 
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Question  1: 
Answer  1: 

Question  2 : 
Answer  2: 

Question  3: 

Answer  3 : 

Question  4: 
Answer  4: 

Question  5: 
Answer  5 : 

Question  6: 
Answer  6: 

Question  7: 
Answer  7: 

Question  8: 

Answer  8: 
Question  9: 

Answer  9: 

Question  10; 
Answer  10: 

Question  11: 

Answer  11: 

Question  12: 

Answer  12: 
Question  13: 


What  should  I  do  if  I  have  a  question  about  this  program 
that  is  not  answered  in  the  documentation? 

Call  Richard  Kraft,  NBS,  Applied  Mathematics  Division  or 
Electron  Devices  Division,   (301)  921-3621. 

Where  is  the  physical  problem  explained? 
In  section  2.1. 

Where  are  the  partial  differential  equations  listed  in  sum- 
mary form? 

Equations   (2.13)   to  (2.19). 

Where  are  the  computational  procedures  outlined? 
In  section  3.1. 

Where  are  the  discretized  equations  derived? 
In  chapters  5  and  6 . 

What  machine  has  the  program  run  on? 
UNIVAC  1108. 

Where  are  the  flow  diagrams? 
In  chapter  8. 

Where  are  the  detailed  explanations  of  the  computational 
steps? 

In  chapter  7. 

Where  are  the  input  and  output  and  a  description  of  how  to 
work  a  simple  problem? 

In  chapter  9. 

Where  is  the  listing  of  the  computational  blocks? 
In  chapter  11. 

Suppose  I  want  to  know  where  to  go  in  the  documentation  to 
learn  about  some  quantity  in  the  listing? 

Go  to  the  glossary  in  chapter  12. 

What  is  the  basis  for  the  program's  validity  and  what  are 
its  limitations? 

See  chapter  10. 

What  about  the  "deck  cards"  and  the  statements  in  the  list- 
ings that  are  not  dociomented  in  VERSION  1? 


Answer  13:        Completely  ignore  them.     Everything  needed  to  understand 

VERSION  1  is  explained  in  this  documentation  and  is  in  the 
computational  blocks   (see  listing,  chapter  11) . 
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2. 


DESCRIPTION  OF  THE  PROBLEM 


2.1.     Physical  and  Mathematical  Description  of  the  Problem 

The  problem  considered  here   [4]  deals  with  the  redistribution  of  dopant 
impurities  found  in  silicon  as  the  silicon  is  oxidized  and  converted 
into  silicon  dioxide.     For  the  sake  of  specificity,  we  consider  the  im- 
purity to  be  boron  and  have  studied  the  following  ideal  situation. 

Suppose  the  entire  half  space  lying  to  the  right  of  the  x-z  plane,  see 
figiire  la,  to  be  occupied  by  silicon  while  the  region  on  the  left-hand 
side  of  this  plane  is  filled  with  oxygen  at  some  elevated  temperature. 
Also,  suppose  that  the  silicon  is  doped  with  boron  and  that  the  density 
distribution  of  this  impurity  is  constant  in  planes  parallel  to  the  x-z 
plane  so  the  distribution  is  essentially  only  a  function  of  the  y- 
dimension. 

The  time  evolution  of  an  initially  uniform  boron  distribution  (for  appro- 
priate conditions  at  the  oxygen-oxide  interface)   is  illustrated  in  figure 
lb.     When  the  oxygen  contacts  the  silicon,  it  begins  to  react  to  form 
silicon  dioxide,  and  a  sharply  defined  moving  boundary  which  separates 
the  pure  silicon  from  that  previously  converted  into  silicon  dioxide  be- 
gins to  move  into  the  pure  silicon.     The  position  of  this  moving  boundary 
is  known  empirically,  see   [1]   and  section  4.2,  and  is  denoted  by 


where  t  is  time. 

Inside  the  silicon  and  oxide,  the  transport  of  boron  is  governed  by  or- 
dinary Fickian  diffusion  but  with  different  diffusion  coefficients  Dj 
and  D2  in  the  oxide  and  silicon,  respectively. 

The  more  complicated  physical  processes  take  place  at  the  moving  oxide- 
silicon  interface.     The  boron  has  a  preference  to  be  dissolved  in  the 
oxide  rather  than  the  silicon.     Hence,  as  the  moving  boundary  moves  into 
the  silicon,  the  boron's  preference  to  be  in  the  oxide  depletes  the  bo- 
ron density  in  the  vicinity  of  the  moving  boundary.     The  preference  of 
the  boron  to  be  in  the  oxide  rather  than  the  silicon  is  quantitatively 
expressed  by  requiring  that  the  moving  boundary  segregation  condition  be 
satisfied. 


y  =  y   (t) ,  for  t  >  0 


(2.1) 


C2(y  (t)  ,t)  = 
o 


m  Ci  (y^(t)  ,t) 


(2.2) 


where  m  is  a  given  constant  (the  segregation  coefficient)  and  eq  (2.2) 
must  be  satisfied  for  all  time  t  >  0.  In  eq  (2.2)  Ci(y,t)  and  C2(y/t) 
represent  the  boron  densities  in  the  oxide  and  silicon,  respectively. 


X 


MOVING  BOUNDARY 


0- 


SiO: 


I  1  ►  Y 

Yo(t) 

Figure  1.     The  problem  geometry. 


In  addition  to  the  condition  eq  (2.2)  holding  at  the  moving  boundary,  it 
is  also  assumed  that  there  is  conservation  of  boron  atoms.  Mathematical- 
ly, this  condition  is  expressed  by  the  requirement  that 


i  y  ^ 


+   (a-^-l)y  Ci 


y=yjt) 


=  DoS  Co 

y  ^ 


y=y^(t) 


(2.3) 


y=y^(t) 


In  order  to  understand  the  origin  of  the  convective  flux  (a~ ■'-1) yQCj ,  it 
is  necessary  first  to  understand  a  second  complicated  physical  effect  oc- 
curring at  the  moving  oxide-silicon  boundary.     When  silicon  is  converted 
into  oxide,  its  volume  expands  by  a  factor  a~^,  where  a      0.45.     In  the 
model  chosen  here   [4] ,  it  is  assumed  that  all  the  oxidation  of  the  sili- 
con takes  place  entirely  at  the  moving  oxide-silicon  interface.     In  this 
case,  the  expansion  tends  to  push  the  previously  formed  oxide  backwards 
at  a  certain  velocity  into  the  oxygen,  creating  another  moving  boundary 
between  the  oxygen  and  oxide,  see  figure  2.     The  velocity  with  which  the 
oxide  is  being  convected  backwards  is  assumed  to  set  up  a  convective 
drift  transport  in  the  oxide's  boron  atoms.     The  magnitude  of  the  veloci- 
ty of  this  convective  flux  of  boron  atoms  is  determined  in  the  following 
manner . 
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BORON  REDISTRIBUTION  DURING  OXIDATION 


TIME 
0 


OXYGEN 


SILICON 


+  1 


+  2 


OXYGEN 


+  3 


OX 


DE 


SILICON 


GOVERNS  ELECTRONIC 
PROPERTIES  OF  SILICON 
DEVICES 

Figure  2.     A  hypothetical  redistribution  process, 


Consider  a  thin  slab  with  width  Ay  of  pure  silicon  located  exactly  at 
the  moving  oxide-silicon  interface  at  some  time        just  before  it  is  con- 
verted into  oxide.     After  some  time  At,  the  moving  boundary  has  trav- 
eled from  one  side  of  this  thin  slab  to  the  other 


Ay  =  y  (t    +  At)  -  y   (t  )  , 

o    o  o  o 


(2.4) 


converting  it  into  a  thin  slab  of  pure  oxide  of  width  a~-^Ay.     The  differ- 
ence between  this  new  larger  width  and  the  original  width  is  assumed  to 
be  the  distance  that  the  previously  formed  oxide   (i.e.,  prior  to  t^)  is 
moved  backward.     Therefore,  this  rate  is 


lim  g-^Ay  -  Ay  _  ^ 
At->0  At  ^  '  dt 


(2.5) 


t=t 


and  the  backward  velocity  of  the  oxide  is  -(a  ^-Dy^  where  y^  = 


dt 


Thus,  the  full  equation  describing  the  transport  of  the  boron  in  the  ox- 
ide is 
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)  Ci(y,t)  =  Di9"Ci  +  (a-l-l)y  (t)8  Ci  , 


(2.6) 


where  the  last  term  in  eq  (2.6)  represents  the  convective  boron  trans- 
port. 

By  introducing  a  coordinate  system  (see  figure  3) ,  moving  with  the  ori- 
gin at  the  oxygen-oxide  interface, 

z  =  y  +   (a-^-l)y  (t) 
o 


or 


(2.7) 


z  =  y  +  z^(t)   -  y^(t)  , 


where 


z  (t)  =  a-^y  (t)  , 
o  o 


(2.8) 


eq  (2.6)   simplifies  to 


3  Ci  =  D.a^Ci  for  0  <  z  <  z   (t)  . 
t-^         -^z-^  —      —  o 


(2.9) 


Moreover,  in  the  z-frame  there  is  no  longer  a  moving  boundary  at  the 
oxygen-oxide  interface  to  contend  with. 

The  moving  boundary  conservation  condition  in  the  z-  and  y-frame  is 

Di3  Ci(z   ,t)   +  z  Ci(z   ,t)   =  D29  C2(y   ,t)   +  y  C2  (y  ,t)    .  (2.10) 
^  z  ^    o  0^0  y        o  0  0 

In  this  model  we  assume  the  transport  of  boron  across  the  oxygen-oxide 
interface  to  be  governed  by 


Di 9  Ci  =  C     (Ci (0,t)   -  C     J  , 
i-~>  out 


(2.11) 


OXYGEN-^ 


OXIDE- 


•SILICON 


Figure  3.     The  moving  z-coordinate  frame. 
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where  Cp  and  Cout         assumed  to  be  known  constants.     Finally,  at  y  =  oo 
we  assume  the  concentration  to  be  stationary  and  equal  to  a  "bulk"  value 


B 


C2(~.t)  =  Cg  . 


(2.11') 


Since  we  intend  to  solve  the  above  equations  on  the  computer,  it  is  manda- 
tory to  truncate  the  infinite  y-domain  by  selecting  some  point  y  =  BDRY 
far  from  the  y-origin  and  requiring  there  that 


C2(BDRY,t)   =  C  . 

B 


(2.12) 


Finally,  at  time  t  =  0  it  is  assumed  that  the  silicon  slab  is  completely 
unoxidized   [y   (0)   =  0]   and  that  the  concentration  has  a  step  function  dis- 
tributxon : 


C2(y,0)  =  C 
C2(y/0)  =  C 


max 


y  <_  NH 

NH  <  y  <  BDRY 


(2.13) 


where  C       ,  NH  >  0  are  parameters, 
max 

In  summary,  we  will  solve  eq  (2.13)  and 


a  Ci  =  Di3  ^Ci  for  0  <  z  <  z  (t)  , 
t  ^         i  2     ^  —     —  o 

)  Co  =  Do9  ^Co  for  y  (t)  <  y  <  BDRY  , 
t^  Y  o       —  — 


Di3  Ci  +  z  Ci 
^  z  ^        o  ^ 


=  DoB  Co  +  y  Co 


z=z  (t) 

o 


y=y  (t) 
o 


(2.14) 
(2.15) 
(2.16) 


C2 (y  /t)  =  mCi (z  ,t)  , 

o  o 


Di3  Ci  -  C   (Ci  -  C  ) 


,z 


=  0  ,  and 


(2.17) 
(2.18) 

(2.19) 


z=0 

C2(BDRY,t)  =  , 

where  the  moving  boundary  yr,(t)  or  z  (t)  =  a~V  (t)  is  a  specified  mono- 
tonic  curve,  see  section  4.2. 
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2.2.     The  Partial  Differential  Equations  in  Their  Equivalent  Integral 
Form 


The  usual  procedure  for  the  numerical  solution  of  partial  differential 
equations  like  eqs   (2.13)  to  (2.19)   comes  through  solving  the  system  of 
algebraic  equations  that  results  from  replacing  the  partial  derivatives 
in  these  equations  by  their  approximating  difference  quotients.  When 
this  standard  procedure  was  tried   (see   [10-13]   for  other  approaches) ,  it 
did  not  lead  to  a  very  accurate  solution.     The  discretization  of  the  fol- 
lowing system  of  equations  equivalent  to  the  partial  differential  equa- 
tions and  which  involves  integrals  in  place  of  partial  derivatives  does 
lead  to  a  system  of  algebraic  equations  that  can  be  solved  readily  to 
yield  an  accurate  solution  of  the  original  system. 

The  system  of  equivalent  equations  is : 


Ci (z,t)dz  =  -Di3^Ci 


+  Di9  Ci 
^  z  ^ 


z=z. 


=  2B 


(2.20) 


for  every  pair  of  points  z^  <  Zg  such  that  0  <_  z^,  Zg  <_  z^{t)  , 


d_ 
dt 


/ 


C2(y,t)dy  =  -029^02 


+  Co 

^  y  ^ 


y=yA 


y=yB 


(2.21) 


for  every  pair  of  points  y^^  <  yg  such  that  yo(t)  f_  y^/  Y-q  —  BDRY, 


Ci(z,t)dz  +  / 


t)dz  +  /        C2(y,t)dy  =  -Di9^Ci 

yo(t) 


+  DoS  Co 

y 


Z=Z7 


,  (2.22) 


y=yB 


C2(y  ,t)  =  m  Ci  (z  ,t)  , 
o  o 


DiS  Ci   -  C   (Ci (0,t)   -  C     ^)  =  0  , 
z  p  out 


C2(BDRY,t)   =  Cg  ,  and 


(2.23) 
(2.24) 
(2.25) 


C2(y/0)  =  C^(y) 


(2.26) 
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where 


C^(y)  = 


'max      y  <_  NH 

:^  NH  <  y  <  BDRY  . 


To  see  that  eq  (2.14)  is  equivalent  to  eq  (2.20)  ,  first  cariry  out  the 
differentiation  on  the  left-hand  side  of  eq  (2.20)  under  the  integral 
sign  and  substitute  for  9^Ci  by  employing  eq  (2.14) .     One  gets 


d_ 
dt 


I      ci  (z,t)dz  =  /  D^a^^c^ 


dz 


(2.27) 


and 


/ 


DiS^  Cldz  =  -Di9^Ci 


z=z. 


z=z. 


(2.28) 


yielding  the  right-hand  side  of  eq  (2.20). 

By  reversing  these  steps  starting  with  the  right-hand  side  of  eq  (2.20), 
one  finds 


/ 


B 


O^Ci   -  Di9  2c,)dz  =  0 
t  ^         ^  z  ^ 


(2.29) 


for  every  0  1.       <  Zg  <_  ZQ(t)  .     Therefore,  the  integral  of  eq  (2.29)  must 
vanish  (it  is  assumed  to  be  continuous)  and  eq  (2.14)  results.     By  using 
trivial  modifications  of  these  same  techniques,  the  remaining  equations 
can  be  shown  to  be  equivalent  to  eqs  (2.13)  to  (2.19). 
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3.     THE  COMPUTATIONAL  PROCEDURE 


3.1.     General  Description  of  the  Computational  Procedure 

The  numerical  approach  to  the  solution  of  the  partial  differential  equa- 
tions in  the  integral  form,  eqs   (2.20)   to   (2.26),  replaces  the  problem  of 
finding  the  continuous  functions  Ci(z,t),  C2(y,t)  with  the  problem  of  de- 
termining the  functions  Ci(z,  tj^) ,  02(7,  tMM)  which  are  only  defined  at 
discrete  points  on  a  finite  difference  grid.     The  time  coordinates  of 
these  grid  points  are  specified  through  the  definition  of  a  sequence  of 
discrete  times  for  MM  =  1 ,  MME  where  MME  is  some  given  positive 

integer  and  the  first  discrete  time,  t^,   is  zero.     The  spatial  coordinates 
of  the  "fixed"  grid  points   (see  sec.  4.1  and  fig.  4  for  further  details) 
are 

z  =   (J  -  l)Az  for  J  =  1,2,...  , 

(3.1) 

y  =   (I  -  l)Ay  for  I  =  1,...,N2  , 

where  N2  is  some  given  integer  index.     The  mesh  widths  Az,  Ay  are  also 
given  quantities  whose  magnitudes  like  those  of  MME  and  N2  need  not  be 
considered  now.     In  addition  to  the  fixed  grid  points  in  eq  (3.1)  ,  there 
are  two  grid  points  located  at  the  moving  boundary   (see  sec.  2.1).  One 
is  at  the  moving  boundary  bordering  the  z-region  and  the  second  is  at 
the  moving  boundary  bordering  the  y-region ,  see  figure  4.     The  space 
and  time  coordinates  of  these  two  moving  boundary  grid  points  are,  re- 
spectively,   (Zq  ( tMM)  '  ^mm)  (yo(tMM)  '^MM)  •  addition  to  seeking  the 
concentrations  at  the  fixed  grid  points  C^ (J,tMM) /  C2(I,tMM)/  we  also  de- 
sire to  find  the  concentrations  at  the  two  moving  boundary  grid  points, 

^1  (^^^^^mm)  '^MM^  '  C2(y^(t     ),t     )    forMM=l,...,  MME. 

o     MM       MM  o     MM  MM 

Because  of  the  initial  condition  eq   (2.26),  these  concentrations  are 
known  at  time  t^  =0. 

Let  us  assume  that  we  have  already  determined  the  concentrations  at  the 

fixed  and  moving  grid  points  up  to  some  arbitrary  time  t^n;  then  we  will 

proceed  to  describe  how  these  same  concentrations  at  the  next  time  t^M  + 
At  =  t^M+i  (where  At  =  tj^^j.  ~  ''^MM^          ^®  found  by  performing  the  opera- 
tions in  the  following  three  steps. 

Step  #1:     Discretize  the  equations. 

By  discretizing  the  equivalent  description  of  the  partial  differential 
equations  in  eqs   (2.20)   to   (2.26)  ,  we  derive  a  system  of  linear  alge- 
braic tridiagonal  equations  that  couple  the  concentrations  at  time  tj^ 
with  those  at  t^M+l*     Specifically,  with  each  of  the  fixed  grid  points 
in  the  z-region,  we  derive  and  associate  an  equation  of  the  form 

ZiCl(J  -  1,TTMMP1)    +  Z2CI  (J,TTMMP1)    +  ZsCKJ  +  1,TTMMP1)    =  Z14  (3.2) 
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Yaxis 


Z  axis 


1/3  and 


Table  1:     FORTRAN  Version  of  Symbols  in  Eqs   (3.2)   to   (3.4)    in  the  Y- 
and  Z-Regions. 


Cl((P  -  l)Az,t^_^^) 

=  CI (P, TTMMPl) 

=  C2(P, TTMMPl) 

t           =  TTMMPl 
MM+1 

o     MM+1  MM+1 

=  CB(1, TTMMPl) 

'^2(yo^Wi^'Wi^ 

E  CB (2, TTMMPl) 

POSITION  OF  MOVING  BOUNDARY 

-H  Ayh- 

3456     789  N2 


I  2 


Y  =  0 


Z  =  0 


MOVING  BOUNDARY 
GRID  POINTS 


-t  1- 


2     3    4  5 


6     7     8     9  10 


Figure  4.  Finite  difference  grids  with  a 
NZO  =  3. 


for  J  =  1,...,  NL  and  where  the  FORTRAN  notation  is  defined  in  table  1. 
Similarly,  with  each  of  the  fixed  grid  points  in  the  y-region,  we  derive 
and  associate  an  equation  of  the  form 

YlC2(I  -   1, TTMMPl)    +  Y2C2 (I, TTMMPl)    +  Y3C2(I  +  1, TTMMPl)    =  Y^  (3.3) 

where  I  =  NR, . . . ,  N2  and  table  1  contains  the  FORTRAN  symbol  definitions. 

The  coefficients  Zi,...,Zn,  Yi  ,  .  .  .  ,Y^.  are  composed  of  the  mesh  width 

quantities  At,  Az,   Ay  and  hence  are  known.     In  addition,   the  quantities 

Zi^,Yi^  contain,  respectively,       (z ,  tj^^jyj)  and  C2(y,tj^^)   and  therefore  these 

equations  couple  the  concentrations  at  the  two  successive  time  steps 

t       and  t       , . 
MM  MM+1 

At  the  moving  boundary  grid  point  in  the  z-region,   there  is  an  addition- 
al equation  called  the  preliminary  equation   (for  reasons  explained  in 
step  #3)  which  expresses  the  value  of       (zq  (tf^+j)  ,  tj^+l)   in  terms  of 
the  concentrations  at  the  neighboring  grid  points,  at  time  t^+i. 
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DM1  CB{1,TTMMP1)    =  DM2  CI (NL,TTMMP1) 


+  DM3  C2 (NR,TTMMP1)    +  DM4  CI (NL  -   1,TTMMP1)  (3.4) 

+  DM5  C2 (NR  +  1,TTMMP1)    +  DM6 

where  the  FORTRAN  symbols  are  defined  in  table  1.     In  this  equation 
DM1,..., DM6  are  composed  of  known  mesh  quantities  and  DM6  contains  con- 
centrations, at  time  tj^^   (at  the  grid  points  coincident  and  adjacent  to 
the  moving  boundary  grid  points) ,  which  are  of  course  assumed  to  be 
known . 

Step  #2:     Put  the  tridiagonal  equations  into  "P,Q"  form. 

By  algebraically  manipulating  eq   (3.2),   it  is  possible  to  put  it  in 
the  f  orm 

C1(J,TTMMP1)  =  P1(J)  C1(J  +  1,TTMMP1)  +  Ql (J)  (3.5) 
for  J  =  1,...,  NL  -  1  and 

CI (NL,TTMMP1)  =  PA(1)  CB(1,TTMMP1)  +  QA(1)  .  (3.6) 
By  similar  manipulations,  eq   (3.3)   can  be  put  into  the  form 

C2(I,TTMMP1)  =  P2(I)  C2(I  -  1,TTMMP1)  +  Q2 (I)  (3.7) 
for  I  =  NR  +  1,...,N2  and 

C2 (NR,TTMMP1)    =  PA(2)    CB(2,TTMMPl)    +  QA(2)    .  (3.8) 

The  "P,  Q"  coefficients  in  these  equations,   i.e.,  PI ( J) ,  Ql ( J) ,  P2(I), 
Q2(I),  PA(L),  QA(L)    for  L  =  1,2  are  algebraic  combinations  of  the  quanti- 
ties Y^,...,  Yi^,   Z^,...,   Zi^  and  are  hence  known. 

Step  #3:     Determination  of  the  concentrations  at  time  TTMMPl  =  t 

MM+1 

The  first  important  use  of  eqs    (3.5)   to   (3.8)   is  that  it  is  possible 
using  these  equations  with  J  =  NL  -  1  in  eq   (3.5)   and  I  =  NR  +  1  in 
eq   (3.7)   to  eliminate  with  the  help  of  eq   (2.23)   all  the  unknown  concen- 
trations on  the  right-hand  side  of  the  preliminary  eq   (3.4).     Thus,  we 
obtain  the  value  of  the  concentration  at  the  moving  boundary  in  the  z- 
region  in  terms  of  known  quantities, 

CB(1, TTMMPl)   =  "Some  function  of   (DMl ,  ,  ,  .  ,DM6 ,   Zi,...,Yi^)"   .  (3.9) 

Since  eq   (3.4)    leads  directly  in  this  way  to  what  will  turn  out  to  be 
the  most  important  concentration  at  time  TTMMPl  in  the  sense  that  with 
this  value  of  CB(1, TTMMPl)   all  the  other  concentrations  are  simply  cal- 
culated, we  call  eq   (3.4)   the  preliminary  equation  for  distinction  pur- 
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poses.  From  this  result  and  the  segregation  boundary  condition  eq 
(2.23)  with  SEG  =  m,  we  can  find  DB (2 ,TTMMP1) . 


CB(2,TTMMP1)    =  SEG   •   CB(l,TTMMPl)    .  (3.10) 

By  using  these  two  concentrations  at  the  moving  boundary  and  back  sub- 
stituting in  eqs    (3.6)   and   (3.8),  we  can  find  C1(NL,TTMMP1)  and 
02 (NR,TTMMPL) . 


Finally,  by  using  these  two  quantities  and  by  back  substituting  in  eq 
(3.5)   starting  with  J  =  ML  -  1  and  by  back  substituting  in  eq  (3.7)  be- 
ginning with  I  =  NR  +  1 ,  we  can  successively  generate  all  the  concentra- 
tions at  the  fixed  grid  points  at  time  t    .  .       TTMMPl . 

MM+i 

By  repeating  steps  one  through  three  starting  at  time  tj^  =  0  when  the 
concentrations  are  known  in  the  y-region   (there  is  no  z-region  since 
the  moving  boundary  starts  at  z  =  0  at  time  t  -  0) ,  all  the  concentra- 
tions at  all  the  grid  points  at  the  discrete  times  t     ,  MM  =  1,...,MME 

MM 

can  be  found. 


3.2.     Short  Guide  to  the  Following  Chapters 


In  this  documentation,  the  derivation  of  the  tridiagonal  linear  eqs 
(3.2),    (3.3),  and   (3.4)   is  carried  out  in  Chapter  5.     The  manipulation 
of  these  equations  into  the  form  eqs   (3.5)   to   (3.9)   is  carried  through 
in  Chapter  5.     The  computational  procedure  is  described  in  detail  in 
Chapter  7,  and  Chapter  8  contains  flow  diagrams  of  the  listing  in  Chap- 
ter 11.     Chapter  9  contains  a  worked  example  and  a  compact  and  precise 
description  of  the  input  and  output  for  the  program.     Chapter  10  de- 
scribes the  validation  and  limitations  of  the  program,  and  Chapters  12 
and  13  contain  a  glossary  for  FORTRAN  symbols  and  references.     The  ap- 
pendices collect  finite  difference  approximations  used  in  Chapter  5  and 
an  "in-house"  plot  documentation. 

The  next  chapter  collects  material  needed  in  Chapter  5  and  gives  a  de- 
tailed description  of  the  finite  difference  grid,  moving  boundary  mo- 
tion, equation  dimensionalization ,  and  plotting  description. 
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4.      DESCRIPTION  OF  AUXILIARY  COMPUTATIONAL  PROCEDURES 


4.1.     A  Detailed  Description  of  the  Finite  Difference  Grid 

In  this  section  we  describe  more  fully  the  finite  difference  grid  al- 
ready introduced  in  section  3.1,   see  eq   (3.1).     We   (see  fig.   4)  begin 
by  describing  how  the  magnitudes  of  the  mesh  widths  Ay,  Az,  At  and  the 
number  of  grid  points,  N2,  on  the  y-axis  are  determined. 

The  mesh  width  Ay  is  calculated  in  DISTRIB  by  dividing  the  program 
input  quantity  BDRY   [see  above  eq   (2.12)]  by  an  integer  input  NYO,  i.e., 

BDRY 

Ay  =    .  (4.1) 

^       NYO  ^  ' 

Consequently,  placing  the  grid  point  index  1  at  y  =  0  and  calling  the 
index  of  the  y-grid  point  at  y  =  BDRY,  N2 ,   implies  that 

N2  =  NYO  +  1  .  (4.2) 

In  order  to  understand  the  motivation  for  the  way  we  define  the  magni- 
tude of  the  mesh  width  Az,  imagine  the  width  of  an  interval  of  size  Ay 
of  silicon  after  it  has  been  converted  into  oxide.     It  swells  to  a  size 
[see  sec.  2,1]  a~-^Ay.     Since  in  DISTRIB  a  is  an  unrestricted  input  and 
is  possibly  very  small,  this  latter  quantity,  a~-'^Ay,  could  be  very  large. 
We  want  to  choose  the  magnitude  of  Az  in  such  a  way  that  there  are 
enough  Az  intervals  in  this  width  to  describe  any  concentration  varia- 
tion in  it.     Therefore,  we  introduce  another  program  input  integer  NZO 
and  use  it  to  define 

Az  =  a" ^ Ay /NZO  (4,3) 

in  such  a  way  that  a~VNZO  %  1  so  that  Az  and  Ay  are  approximately  the 
same  magnitude.     Another  very  important  advantage  of  this  choice  of  Az 
is  that  it  implies  that  when  the  moving  boundary  sweeps  through  a  mesh 
width  of  magnitude  Ay  in  the  y-region  NZO  mesh  widths  of  magnitude  Az 
are  swept  out  in  the  z-region. 

In  order  to  explain  the  reasoning  behind  the  definition  of  At,  we  re- 
mark that  the  computer  program's  numerical  stability  seems  to  be  af- 
fected if  the  moving  boimdary  moves  too  rapidly  between  successive  grid 
points.     Therefore,  At  is  defined  in  such  a  way  that  the  moving  boundary 
does  not  travel  more  than  a  fraction  WFRC   (real  program  input)  of  a  Az 
mesh  width  between  the  times  t  =  0  and  t  =  At.     Thus,  we  solve  the 
equation 

WFRC  •  Az  =  z   (At)   -  z   (0)  (4.4) 
o  o 

for  At  and  then  define  At  as  this  solution.     This  is  how  At  for  the 
first  iteration  (see  sec.   3.1)  when  t^  =  0  and  Zq(0)  =  0  is  chosen.  In 
computational  blocks  #16  and  #20   (see  sec.   7.1),   it  is  explained  how  At 
is  increased  during  later  iterations  of  the  main  loop  78  in  DISTRIB. 
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time=  t^iM  +  i  ^ 


time  =  tMM  • 


NLMI  NL 
SM(I)AZ 


MOVING  BOUNDARY 
GRID  POINTS 


NLMI  NL 

Figure  5 :     Motion  of  the  moving  boundary  in  the  z-region . 
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Figure  6.     Motion  of  the  moving  boundary  in  the  y-region. 


Before  leaving  the  description  of  the  finite  difference  grid,  it  is  im- 
portant and  convenient  to  introduce  special  names  for  certain  indices 
of  special  grid  points.     Figure  5  shows  the  motion  of  the  moving  boun- 
dary starting  at  time  t|yjj,,j  and  ending  at  time  t]\^+i  in  the  z-region.  We 
denote  the  index  of  the  z-grid  point  which  is  at  or  immediately  to  the 
left  of  the  moving  boundary  at  time  by  NL  and  the  next  grid  point 

to  its  left  by  NLMI   (=  NL  -  1) ,     Similarly   (see  fig.   6) ,  the  index  of 
the  first  grid  point  in  the  y-region  greater  than  the  position  of  the 
moving  boundary  at  time  tj^^  is  denoted  by  NR,  and  the  index  of  the  next 
point  to  its  right  is  NKPl   (=  NR  +  1) .     Of  course,   the  values  of  NL  and 

NR  change  with  t,„  . 

MM 

4.2.     The  Description  of  the  Moving  Boundary 

In  DISTRIB  the  user  can  choose  between  two  different  formulas  for  the 
moving  boundary.     These  are  the  two  motions  related  to  the  wet  and  dry 
oxidations  in   [1].     In  both  cases  there  are  three  parameters  A,  B,  and 
T   (each  with  values  that  are  temperature  dependent)   in  the  equations 
that  describe  the  motion  of  the  moving  boundary. 

Case  I:     Wet  oxidation,  t  =  0. 
The  equation  is   [see  1,  eq  (13)] 


A  r  4tBll/^ 


-  1  for  t  >  0 


(4.5) 
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The  equation  expressing  t  in  terms  of  is 

t  =    (z^^  +  Az^)/B  .  (4.6) 

Table  1  in   [1]  gives  the  values  of  A  and  B  for  different  temperatures. 
Case  II:     Dry  oxidation,   t  >  0. 
For  times  t  >  x  [1] 

■/2 

=  -  I  1  +  ipr-\        ~  ^         ^  >.      -  ('i-'^) 

When  t  <  T ,  we  use  [1] 


A  r  4tBl^' 


where 


ZCRIT  ^ 

z     =  —         •  t  ,  4.8 

O  T 


jl/2 

ZCRIT  =  CNSTX  =  ^  11  +  -  1   .  (4.9) 


A  r  4tBT 


CNSTX  is  a  duplicate  name  for  ZCRIT.  Table  II  in  [1]  gives  the  values 
of  A,  B,  and  T  for  different  temperatures. 

The  FORTRAN  symbols  for  the  constants  in  the  previous  equations  are 
given  in  table  2. 


Table  2:     FORTRAN  Version  of  Symbols  Used  in  the  Moving  Boundairy  Equa- 
tion and  Finite  Difference  Grid. 

A  =  CNSTA 

B  =  CNSTB 

. 5A  =  CNA 

T  =  CNTAU  =  TLIN 

ZCRIT  =  CNSTX 

Az  =  DZ 

Ay  =  DY 

At  =  DT 

a  =  ALPHA 
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4.3.  Dimensionalization 


Although  the  physical  input  in  DISTRIB  is  required  to  be  in  units  of 
micrometers  and  hours,   it  is  computationally  more  efficient  to  solve 
the  partial  differential  equations  with  the  physical  quantities  ex- 
pressed in  other  units.     For  the  kinds  of  problems  of  interest  in  im- 
purity redistribution,  a  convenient  basic  length  is  BLTH  =  0,02  ym. 
Given  the  basic  length  unit,  a  basic  time  is  frequently  introduced  in 
diffusion  theory  by  defining 

2 

BTM  =  BLTH_  j^Q^j^g  (4.10) 
D 

o 

where  Dq  is  some  diffusion  coefficient  in  the  problem.     In  our  program, 
Dq  is  chosen  to  be  D2,  the  diffusion  coefficient  in  silicon.     When  these 
new  units  are  used  to  dimensionalize  the  diffusion  equation  in  the  sili- 
con, one  finds  the  new  equation  has  a  diffusion  constant  of  unity.  That 
is ,  making  a  change  of  independent  variables 

/  _    t        /  _      y  ,     /  ^  z 

^     -  BTM'   ^     ~  BLTH'  ^  BLTH 

transforms  eq   (2.14)  into 

8   ,  C2  =  C 


2  -  °   /  ^2 

t  y 


and  eq   (2.15)  into 


\'  ^1  =      'I'  ^1  • 

The  rest  of  the  dimensionalizations  carried  out  in  block  #4  of  section 
7.1  are  standard  conversions  from  micrometers  and  hours  to  the  new 
units  BLTH  and  BTM. 

4.4.     The  Description  of  the  Way  Arrays  CI  and  C2  Are  Plotted 

The  plotting  is  really  a  peripheral  item  and  not  described  in  as  much 
detail  as  substantive  parts  of  the  program.     All  the  quantities  un- 
defined are  defined  or  referenced  in  the  glossary  and  their  definition 
in  the  glossary  makes  their  function  clear.     First,  the  way  CI  is 
plotted  is  described.     Let  NC  be  the  integer  part  of   (NL  -  1)/NM0D. 
When  NC  >  0,  we  plot  the  concentrations  CI  and  z-coordinate   (in  units 
of  ZOPl)   at  every  grid  point  K  =  V  •  NMOD  +  1  for  V  =  1,... ,  NC.  When 
NC  =  0,  we  plot  the  concentrations  at  z  =  0  and  at  the  moving  boundary 
Zq.     The  quantities  XlS,  YlS  determine  the  largest  abscissa  and  ordi- 
nate on  the  graph  of  CI  versus  z. 
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The  graph  of  array  C2  is  plotted  in  the  following  way.     All  abscissa 
are  in  units  of  UN2   (which  is  essentially  SCALE2,  since  y  is  still  in 
units  of  BLTH) .     The  first  C2  concentration  plotted  is  CB(2)   at  the 
moving  boundary.     The  next  concentration  is  C2 (NR)  with  abscissa  S(2)DY 
Next,  we  plot  ND  -  2  concentrations  with  the  indices  KS  =  NR  +   (V  -  2) 
NMODl  for  V  =  3,...,  ND.     Next,  we  pass  these  grid  points  through  a 
filter  which  eliminates  those  with  an  abscissa  greater  than  X2S/2  = 
X2SCT.     The  quantities  X2S  and  Y2S  determine  the  largest  abscissa  and 
ordinate  that  appear  on  the  C2  versus  y  plot.     Note:     The  specific 
scales  used  as  the  ordinates  and  abscissa  for  these  plots  were  chosen 
in  order  to  make  the  resultant  graphs  identical  in  size  to  those  ob- 
tained from  the  equipment  which  experimentally  determined  the  concen- 
tration data.     This  also  accounts  for  the  elimination  of  points  beyond 
X2S/2.     Both  conditions  may  be  changed  to  accommodate  the  user. 

The  smallest  abscissa  and  ordinate  are  determined  in  both  graphs  by 
(0,0) . 

The  plot  subroutine  is  an  NBS  "house  program"  that  is  documented  in  Ap- 
pendix 2. 

This  plotting  procedure  has  worked  well  in  all  cases,  the  only  adjust- 
ments needing  to  be  made  being  in  NMOD  and  NMODl   (see  table  20  in 
sec.  9.22). 
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5.      DERIVATION  OF  THE  DISCRETE  ALGEBRAIC  EQUATIONS 


5.1.     The  General  Prescription  for  Obtaining  the  Tridiagonal  Equation 
Associated  with  a  Fixed  Grid  Point 


This  section,  which  should  be  read  casually,  gives  a  rough  description 
of  the  precise  derivations  in  the  following  sections. 

Each  particular  tridiagonal  equation  in  eqs    (3.2),    (3.3)    [and  (3.4) 
also]  associated  with  some  particular  grid  point  is  derived  by  discre- 
tizing  some  integral  form  of  the  conservation  of  number  eqs   (2.20)  to 
(2.22)   over  the  "flux  cell"   (to  be  described  below)   that  contains  that 
particular  grid  point.     The  flux  cells  are  intervals  that  arise  by  par- 
titioning the  z-  and  y-regions  in  a  certain  way.     That  is,  both  the  z- 
region  and  the  y-region  are  broken  up  into  nonover lapping  intervals 
called  flux  cells  or  concentration  cells  in  such  a  way  that  the  inter- 
vals completely  cover  the  regions   [0,Zo(t)]   and   [yo(t),BDRY]   at  every 
time  t;  and  each  grid  point  is  contained  in  at  most  one  cell  except  for 
the  grid  points  with  indices  NL  and  NR.     Each  of  these  grid  points  is 
contained  in  both  the  irregular  sized  cell  with  which  it  is  associated 
and  the  "moving  boundary  cell";   these  statements  will  be  explained  more 
completely  later.     The  two  boundaries  of  each  of  the  concentration  cells 
are  called  the  flux  boundaries  of  the  cell.     We  interpret  the  concentra- 
tion, C1(J,TTMMP1)   =  C1((J  -  l)Az,ti4M+i),   at  some  grid  point  J  to  repre- 
sent the  average  number  density  in  the  cell  at  time  t^^+j. 

The  flux  boundaries  for  the  cells  in  the  y-region  are  shown  in  figure  7 
by  the  vertical  lines  located  halfway  between  successive  grid  points. 
Notice  that  all  the  cells  in  the  y-region  that  contain  fixed  grid  points 
(see  sec.   3.1)   are  of  width  Ay  except  the  irregular  one  of  width  0.5  Ay 
that  contains  and  leads  to  the  tridiagonal  algebraic  equation  associated 
with  the  grid  point  with  index  NR  at   (NR  -  l)Ay.     This  cell  has  one  flux 
boundary  at   (NR  -  l)Ay  and  the  other  is  at   (NR  -  1/2) Ay.     Similar  re- 
marks apply  to  figure  8. 

The  significance  of  the  flux  cell  containing  some  grid  point,  with  in- 
dex I  in  the  y-region   (see  fig.   7)   is  that  it  leads  to  a  linear  alge- 
braic equation  containing  the  concentration  C2(I,TTMMP1)   as  an  unknown 
when  the  conservation  eq   (2.21)   defined  on  this  cell  is  discretized. 


NR  NRPI 


I-l 


lime  =  t^M  +  i 
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N2 
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1  N2 

f  1  ■  1  ■  1 

■     1     ■     1     ■     1  ■ 
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t 
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Figure  7:     Flux  cells  and  flux  boundaries  in  the  y-region  for  fixed 
grid  points. 
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time  =  tMM+i 


Figure  8.     Flux  cells  for  fixed  grid  points  in  the  z-region  when  NL  > 


In  particular,   let  yj^  and  yg  in  eq   (2.21)   be  at   (I  -  3/2)  Ay  and 

(I  -  1/2) Ay,  respectively.     Make  now  the  following  approximations  in 

eq   (2.21)    (with  t  =  t^^.  )  : 

MM+1 


d 

dt 


f 


C2(y,t)dy  ^ 


C2(y,t  )]dy/At 
MM 


^2(y'Wl^^^  '^^  C2((I 


/ 


C2(y't^)dy      C2((I  -  l)Ay,t^)Ay  , 


Do8  Co 

y 


D2[C2((I  -  l)Ay,t^^^  )   -  C2((I  -  2)Ay,t^^^  )]    ,  and 


y=yA 


C2 


Z  D2[C2(lAy,  t^^^)   -  C2((I  -  DAy,   t^^^^  )  ]  . 


MM+1 


y=yB 


Then  eq   (2.21)   becomes,  up  to  terms  of  higher  order,   the  linear  tridi- 
agonal  equation  in  eq   (3.4)   associated  with  the  grid  point  with  index  I 


C2((I  -  2)Ay,t^^^)   -  |2  -H  C2((I  -  l)Ay,t^^i) 

+  Co(lAy,t         )   =  -  C2((I  -  l)Ay,t  ) 

MM+1  n^A-i-     ^  MM 


D2At 
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The  important  fact  this  example  illustrates  is  that,  by  evaluating 
the  fluxes  in  eq  (2.21)   at  the  flux  boundaries  of  the  flux  cell  and 
evaluating  the  integrals  by  a  quadrature  approximation  with  a  "node" 
located  at  the  grid  point  in  the  flux  cell,  one  ends  up  with  a  tridi- 
agonal   (i.e.,  an  equation  involving  the  concentrations  at  three  con- 
secutive grid  points)  equation  for  the  unknown  concentration  at  time 
''^MM+1        that  grid  point.     Thus,  with  the  help  of  these  flux  cells  one 
can  in  a  natural  way  associate  with  each  unknown  at  a  fixed  grid  point 
an  algebraic  equation. 

The  flux  cells  of  the  points  in  the  z-region  depend  on  the  values  of  NL. 
There  are  three  cases  depending  on  whether  NL  =  1,  NL  =  2,  or  NL  >  2. 
The  cell  configuration  for  the  case  NL  >  2  can  be  visualized  from  figure 
8.     In  the  case  NL  =  2   (see  fig.   9) ,  there  are  just  two  irregular  cells 
each  of  width  0.5  Az.     The  flux  cell  associated  with  the  grid  point  with 
index  NL  has  its  right  side  flux  boundary  located  coincident  with  the 
grid  point.     The  flux  cell  associated  with  the  first  z-grid  point  with 
index  1  has  its  left-hand  flux  boundary  coincident  with  the  first  grid 
point.     These  facts  also  hold  when  NL  >  1,  but  in  this  case  there  are 
regular  grid  cells,   i.e.,  with  width  Az  in  between  the  irregular  ones 
at  the  outer  boundaries.     The  case  NL  =  1  is  special  because  there  are 
no  grid  points  in  the  z-region  between  the  moving  and  fixed  boundaries. 

Finally,  we  point  out  that  for  all  values  of  NL  the  conservation  formula 
eq   (2.22)    (conservation  across  the  moving  boundary)   is  discretized  to 
derive  the  preliminary  eq  (3. 4),  and  in  this  derivation  the  flux  cell  is 
as  shown  in  figure  10  with  the  flux  boundaries  at  z^  =   (NL  -  l)Az  and 
yg  =   (NR  -  l)Ay.     The  trapezoidal  rule  is  used  to  evaluate  the  integrals 
in  eq   (2.22)  with  the  nodes  at  NL,  NR  and  the  moving  boiindary  grid 
points  at  S(l)Az  +   (NL  -  l)Az  and   (NR  -  l)Ay  -  S(2)Ay,   see  figures  5 
and  6. 

These  different  cases,  which  arise  because  we  have  flux  cells  with  dif- 
ferent widths  and  with  different  orientation  of  grid  points  relative  to 
the  flux  boundaries  for  different  values  of  NL,  explain  the  slight  vari- 
ations, classified  in  tables  3  and  15,  in  the  way  the  tridiagonal  alge- 
braic equations  and  their  "P,Q"  coefficients  (see  sec.  3.1)  are  derived 
in  the  following  sections. 


NL=2 


1^  \  f 


time  =  +1 


time  =  t 


MM 


Figure  9.     Flux  cells  for  fixed  grid  points  in  the  z-region  when  NL  =  2. 
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Figure  10.     Schematic  description  of  the  flux  cell  associated  with  the 

discretization  of  the  conservation  equation  across  the  moving 
boundary . 


Table  3:     Guide  to  the  Sections  Where  the  Tridiagonal  Algebraic  Equations 
Associated  With  the  Grid  Points  are  Described  and  Derived. 


Case  I:       The  moving  boundary  has  not  reached  the  second  z-grid  point  or 
NL  <  2. 

(a)  The  index  of  the  z-grid  point  is  1   (see  sec.  5.9) 

(b)  The  index  of  the  y-grid  point  is  greater  than  NR  (see  sec. 
5.3) 

(c)  The  index  of  the  y-grid  point  is  NR  (see  sec.  5.5) 

(d)  The  grid  point  is  at  the  moving  boundary   (see  sec.  5.7) 

Case  II:     The  moving  boundary  has  reached  the  second  z-grid  point  or 
NL  ^  2. 

(a)  The  index  of  the  z-grid  point  is  1   (see  sec.  5.8) 

(b)  The  index  of  the  z-grid  point  ranges  from  2  to  NL  -  1  (see 
sec.  5.2) 

(c)  The  index  of  the  z-grid  point  is  NL   (see  sec.  5.4) 

(d)  The  grid  points  are  at  the  moving  boundary   (see  sec.  5.6) 

(e)  The  index  of  the  y-grid  point  is  greater  than  NR  (see  sec. 
5.3) 

(f)  The  index  of  the  y-grid  point  is  NR  (see  sec.  5.5) 


5.2.     The  Index  of  the  Z-Grid  Point  Ranges  from  2  to  NL  -  1  When  NL  >  2. 

In  this  case  we  suppose  the  index,  J  +  1,  is  such  that  2<J+1<NL-1, 
and  begin  with  the  formula,  see  eq  (2.20) 
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t)dz  =  -  Di 3  Ci 
z  ^ 


(5.1) 


For  the  z-grid  point,  J  +  1,  we  assume  zp^  and  Zg  "the  flux  boundaries" 
to  be  at   (J  -  1/2) Az,    (J  +  1/2) Az,  respectively.     We  speak  of  these  grid 
points  as  regular  grid  points  because  the  distance  betv/een  the  cell  flux 
boundaries  is  Az. 

The  left-hand  side  of  eq  (5.1)   is  discretized  by  employing  the  backward 
time  difference  formula   (A7)   in  Appendix  1  between  the  times  t  = 
and  t  =  tiyu4  +  At.     Here,  At  is  some  small  positive  time  increment  that 
sometimes  depends  on  the  index  MM.     After  the  derivative  has  been  re- 
placed by  the  backward  difference  quotient,  the  integrals  are  approxi- 
mated by  the  midpoint  quadrature  rule   (A2) . 

Finally,  after  further  algebraic  manipulation  we  have 

•  z 


^  Ci(z,t)dz  =  II  [Ci(JAz,t^  +  At) 


z. 


-  Ci(JAz,tj^)]   +  0(Az2  +  At2)/At  .  (5.2) 

Since  all  the  expressions  in  eq  (5.1)  are  assumed  to  be  evaluated  at 
t  =  t]yiM  +  At,  the  approximation  of  the  derivatives  on  the  right-hand 
side  of  eq  (5.1)  gives 


9  Ci 
z 


and 


9  Ci 
z  ^ 


=   [Ci  (JAz,t      +  At)   -  Ci((J  -  l)Az,t,^^  +  At)]/Az  +  0(Az)  (5.3) 
J-  MM  MM 


z=ZA 


=    [Ci((J  +  l)Az,t,,,  +  At)   -  Ci(JAz,t^„  +  At)]/Az  +  0(Az)  .(5.4) 
MM  MM 


Z=Zt 


By  substituting  eqs   (5.2)   to   (5.4)   into  eq   (5.1)   and  further  rearranging, 
we  obtain  up  to  higher  order  terms  the  discrete  version  of  eq  (5.1)  in 
tridiagonal  form, 


DiCi((J  -  l)Az,t^  +  At)   -    (2   •  D    +  Az2/At)   Ci(JAz,tj^  +  At) 

+  DiCi  (  (J  +  l)Az,t^  +  At)   =  ^-^^^  Ci  (JAz,t^) 
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After  introducing  the  FORTRAN  quantities  in  table  4,   eq   (5.5)  becomes 
DFl    •   C1(J,TTMMP1)    -    (2    •   DFl  +   DZ**2/DT)    •  CI ( J  +  1,  TTMMPl) 

(5.6) 

- (DZ) **2 

+  DFl    •   C1(J  +  2), TTMMPl)    •   CI ( J  +  1,TTMM) 

Table  4:     FORTRAN  Version  of  Quantities  in  Eq  (5.5). 
Di   =  DFl 

Cl(    ,    )   =  Cl(   ,  ) 
=  TTMM 

MM 

t       +  At  =  TTMMPl 
MM 

At  =  DT 

Az  =  DZ 

Ci(JAz,t.„)   =  CI  (J  +  1,TTMM) 
MM 

Ci  (JAz,   t       +  At)   =  C1(J  +  1,  TTMMPl) 

MM 


In  the  FORTFIAN  program,  the  following  intermediate  quantities  in  table  5 
have  been  introduced  for  convenience  and  computational  efficiency. 


Table  5:     Intermediate  Quantities  Used  in  Eq  (5.7). 


DFl  =  A2L 

-(2    •   DFL  +    (RAl)"^)    =  A2M 

DFl  ^  A2R 

- (RA1)~^   =  AIM 

AF(J  +  1)    =  AIM   •   C1(J  +  1,TTMM) 


In  terms  of  these  intermediate  quantities,  eq  (5.6)   can  be  written  in 
the  compact  tridiagonal  FORTRAN  form 

A2L   •   CI (J, TTMMPl)    +  A2M   •   Cl ( J  +  1, TTMMPl)    +  A2R  •   CI ( J  +  2, TTMMPl) 

(5.7) 

=  AF  (J  +  1)  . 

As  has  been  explained  in  more  detail  in  section  3.1,  the  quantities 
A2L,...,  AF(J)  are  used  in  computing  the  values  of  the  arrays  PI,  Ql, 
and  Dl.     Specifically,   the  quantities  A2L, . . . ,  A2R,  AIM  are  computed  in 
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subroutine  INTER.     The  subroutine  INTER  is  called  by  subroutine  TRIDNL 
when  this  latter  subroutine  is  calculating  the  arrays  PI,  Ql.     The  array 
AF(J)   is  computed  in  block  #25  of  the  main  program,  see  section  7.1. 

SUMMARY 

Schematic  summary  of  the  derivation  of  the  tridiagonal  FORTRAN  equation 
for  z-grid  points  with  indices  between  1  and  NL: 

(5.1)  ->  (5.5)  (5.6)  (5.7)  . 

Summary  notes: 

(1)  Quantities  in  table  5  except  for  AF  are  computed  in  subroutine 
INTER. 

(2)  The  array  AF (J)   is  computed  in  block  #25  of  the  main  program. 


5.3.     The  Index  of  the  Y-Grid  Point  is  Greater  than  NR 


The  derivation  of  the  tridiagonal  FORTRAN  form  of  the  algebraic  equa- 
tions associated  with  these  grid  points  is  almost  entirely  analogous 
to  the  derivation  of  the  equations  in  siibsection  5.2  once  the  substitu- 
tions z  for  y,         for  D^,         for  C^,  etc.,  are  made.     However,  there  is 

one  exception;  in  the  y- region  we  have  divided  the  equivalent  of  eq 
(5.5)  by  D2  so  that  the  equivalent  of  that  equation  is 

C2((J  -  l)Ay,t  _  +  At) 


MM 


2  + 


Ay2_] 
D2AtJ 


D2At 

+  C2((J  +  l)Ay,t 


C2(JAy,t„„  +  At) 


MM 


MM 


+  At)  = 


5f^C2(JAy,t^) 


(5.8) 


Table  6  is  the  equivalent  of  table  4. 


Table  6:     FORTRAN  Version  of  Quantities  in  Eq  (5.8). 

D2  =  DF2 
Ay  =  DY 

C2(-,-)   =  C2(-,-) 

Az  =  DZ 

C2(JAy,t,,)   =  C2(J  +  1,TTMM) 
MM 

C2(JAy,t       +  At)   =  C2(J  +  1,TTMMP1) 
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When  the  quantities  in  table  6  are  introduced  into  eq   (5.8),   it  becomes 

[DY**2  "I 
2  +  ^ff~r~^  \  C2(J  +  1,TTMMP1)    +  C2(J  +  2,TTMMP1) 


(5.9) 
DY**2 

=   C2  (J  +  1,TTMM)  . 

DF2    •  DT 

This  form  of  the  equation  is  further  simplified  by  introducing  the  inter- 
mediate quantities  in  table  7. 


Table  7:     Intermediate  Quantities  Used  in  Eq  (5.10). 
DT 

 ^  =  RA2 

DY2 

1  =  B2L 

-(2  +  RA2~VdF2)    =  B2M 
1  =  B2R 

-RA2"VdF2  =  BIM 

BF(J  +  1)    =  BIM   •    C2 ( J  +  1,TTMMP1) 


In  terms  of  the  intermediate  quantities  in  table  7,   eq   (5.9)  assumes 
the  form 

B2L  C2(J,   TTJ4MP1)    +  B2M  C2  (J  +  1,   TTMMPl)    +  B2R  C2  (J  +  2,  TTMMPl) 

(5.10) 

=  BF(J  +  1)  . 

The  quantities  RA2 ,   B2L,...,B1M  are  computed  in  subroutine  INTER.  This 
subroutine  is  called  by  subroutine  TRIDNL  and  in  this  routine  the  quan- 
tities RA2,...,B1M  are  used  to  compute  arrays  P2 ,  D2.     The  quantity 
BE (J)   is  computed  in  block  #25  of  the  main  program. 

SUMMARY 

Schematic  siammary  of  the  derivation  of  the  tridiagonal  FORTRAN  equations 
for  y-grid  points  with  indices  greater  than  NR: 

(5.8)  — >  (5.9)  — >  (5.10)  . 

Sammary  notes: 

(1)  Quantities  in  table  7  except  for  BE  are  calculated  in  sub- 
routine INTER. 

(2)  The  array  BE  (J)   is  compuced  in  block  #26  of  the  main  program. 
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5.4.     The  Index  of  the  Z-Grid  Point  is  NL  and  NL  >  1 


In  this  subsection  we  derive  the  FORTRAN  tridiagonal  form  of  the  equa- 
tion associated  with  the  z-grid  point  with  index  NL  which  is  located  at 
z  =   (NL  -  l)Az. 

The  derivation  begins  with  the  integral  conservation  equation 


d 

dt 


/ 


-—  /         Ci(z,t)dz  - 


A 


-Di9  Ci 
^   z  ^ 


Z=Z7 


Di9  Ci 


(5.11) 


Z=Zt 


In  this  equation,   Zg  =    (NL  -  l)Az,   Zji^  =   (NL  -  3/2)  Az  are  the  "flux 
boundaries";  notice  the  flux  boundary  at  Zg  coincides  with  the  location 
of  the  grid  point  associated  with  the  cell   (see  fig.   9).     Also,  the 
cell  width,  the  distance  between  the  flux  boundaries,  is  Az/2 ,  and 
therefore  this  grid  point  is  called  an  irregular  grid  point,  see  sec- 
tion 5.1. 

By  using  the  backward  difference  quotient  approximation  in  eq  (A7)  to 
approximate  the  time  derivative  on  the  left-hand  side  of  eq  (5.11)  be- 
tween the  times  tjyuyi  and  ty^j^  +  At  and  the  right  end  point  quadrature  ap- 
proximation in  eq  (A3)   to  approximate  the  resulting  integrals,  the 
left-hand  side  of  eq  (5.11)  becomes  equivalent  to 


d_ 
dt 


/ 


(NL-1) Az 


Ci(z,t)dz  =  ff^  ]ci[(NL  -  DAz,  +  At] 


(NL-3/2)  Az 


(5.12) 


-  Ci[(NL  -  l)Az,  t^]     +  0(At  +  AzVAt)|  . 

The  spatial  derivative  on  the  right-hand  side  of  eq  (5.11)  at  z  =  z^  is 
approximated  by  the  central  difference  approximation 


'd  Ci 


z=z 


|^|ci[(NL  -DAz,   t^  +  At]   -  Ci[(NL  -  2)Az,   t^  +  At]| 


A 


(5.13) 


+  0(Az)  . 


To  approximate  the  derivative  at  z  =  zg  on  the  right-hand  side  of  eq 
(5.11),  we  use  the  three  point  derivative  approximation  given  in  eq 
(A9) .     Setting  in  that  formula 
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XI  =  z   (t      +  At) 
o  MM 

X  =    (NL  -  l)Az  =  z 
X2  =    (NL  -  2)Az 
h  -  Az 


B 


we  obtain  the  sought-after  result 

1  ( 


9  Ci 
z  ^ 


S  (1) 

Al  )-  1  +  S(l)   ^1^^^^  -  2)^-'  +  At] 


z=z 


B 


+    [(S(l)   -  1)/S{1)]   Ci[NL  -  DAz,   t^^^  +  At] 

MM 


(5.14) 


where 


+  — 7T-^ — rr^  TTTT  Ci  [z    (t       +  At)  ,    t       +  At]  }  +  0(Az2) 

S(l)     [1  +  S(l)]      ^     o     MM  MM  I 


S(l)   =    [z    (t^^^  +  At)    -    (NL  -  l)Az]/Az 
o  MM 


(5.15) 


When  eqs  (5.12),  (5.13),  and  (5.14)  have  been  substituted  into  eq  (5.11) 
and  the  expression  simplified,   it  becomes  up  to  terms  of  higher  order 

DiS(l) 

-  1  +  S(l)   ^1^^^^  -  2)^^'  +  At] 


r^i  •  s(i) 


R(l) 


Ci  [  (NL  -   l)Az,   t       +  At] 
MM 


(5.16) 


(1  +  S(l))    ^if^o^^MM  ^MM  + 


S(l)   Ci  [(NL  -  DAz,  t^^] 

MM 


R(l) 


where 


R(l)  = 


2At 

Az2 


(5.16') 


After  introducing  the  FORTRAN  notation  in  table  8  into  eq  (5.15),  it 
becomes 
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DFl  DFl  DFl 

C1(NL  -   1,   TTMMPl)    +      ,\    +      ,\    CI  (NL,TTMMP1) 


1  +  S(l)  R(l)  S(l) 

(5.17) 

DFl  CB(1)    -  CKNL,  TTMM) 


S{1)     (1  +  S(l))  R(l) 

Table  8:     FORTRAN  Version  of  Quantities  in  Eq  (5.16). 

Di   =  DFl 

t       =  TTMM 
MM 

t.  .    +  At  =  TTMMPl 
MM 

At  =  DT 
Az  =  DZ 

Z    •  DT 


R(l)  - 


DZ**2 


Ci((NL  -   2)Az,    t       +  At)    =  CKNL  -   1,  TTMMPl) 
^  MM 

Ci((NL  -   l)Az,    t       +  At)    -  CKNL,  TTMMPl) 
■L  MM 

Ci((NL  -   l)Az,    t^)    =  CKNL,  TTMM) 

Cl(z    (t       +  At),   t       +  At)   =  CB(1) 

■L     o  MM  MM 


Equation   (5.17)  will  be  used  to  obtain  the  "P,  Q"  quantities  PA(1), 
QA(1)   for  z-grid  point  NL,   i.e.,  quantities  such  that 

CKNL,   TTMMPl)    =  PA(1)    •   CB(1)    +  QA(1)    ,  (5.18) 

in  section  6.6. 

SUMMARY 

Schematic  summary  of  the  derivation  of  the  tridiagonal  FORTRAN  equation 
for  the  z-grid  point  with  index  NL  when  NL  >_  2 , 

(5.11)  — ^   (5.16)  — >  (5.17)  . 

Summary  note : 

Equation  (5.17)  is  used  in  section  6.6  to  obtain  quantities  PA(1) , 
QA(1)  . 

5.5.     The  Index  of  the  Y-Grid  Point  is  NR 

The  derivation  of  the  tridiagonal  FORTRAN  formula  for  this  grid  point 
is  analogous  to  the  derivation  in  section  5.4  for  the  irregular  z-grid 
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point  at  NL  so  the  user  should  refer  to  that  section  if  more  details 
concerning  the  following  derivation  are  desired. 


Again,  we  start  with  the  integral  conservation  equation. 


•4 


B 


C2(y,t)dy  =  -D23  C2 

y 


+  Do9  Co 

y 


y=yA 


(5.19) 


y^yB 


where  y    =   {NR-l)Ay  and  y     =   (NR  -  i)Ay. 

By  using  the  same  approximation  procedure  that  led  from  eq   (5.11)   to  eq 
(5.12)    (except  for  using  the  left  end  point  quadrature  approximation), 
we  find 


dt  /     ^  C2(y^t)dy  =  1^  jC2[(NR  -  l)Ay,  t^  .  At] 


-  C2  [  (NR  -  1)  Ay,   t  ] 
^  MM 


(5.20) 


The  spatial  derivative  on  the  right-hand  side  of  eq   (5.19)   at  y  =  yB  is 
approximated  by  the  central  difference  approximation   (see  eq  (A8) , 
with  h  =  Ay/2) 


3  Co 
Y  ^ 


Ay 


C2(NRAy,   t^^  +  At)   -  C2  [  (NR  -  l)Ay,  t^^  +  At])-+  0(Ay) 


MM 


y=yB 


;5.2i) 


The  spatial  derivative  at  y  =  y^  on  the  right-hand  side  of  eq  (5.19)  is 
approximated  by  using  the  three  point  formula  eq  (A9)  with: 


We  find 

8  <Z. 


XI  =  NRAy 
X  =    (NR  -  l)Ay 

X2  =  y  (t^„  +  At)  . 
o  MM 


T-  It-^-^^W  C2((NR)Ay,  t  +  At) 
AyJl  +  S(2)     ^  ^  MM 


y=yA 


.         -(2)'^^   C2((NR-  DAy,   t^  .  At) 


(5.22) 


S(2)(l  ^  S(2))    ^^(^o^^MM  ^  ^MM  ^  °  ^'^^  ' 
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where 


S(2)   =  -[yo(t^  +  At)   -    (NR  -  l)Ay]/Ay  .  (5.23) 


When  eqs  (5.20),  (5.21),  and  (5.22)  have  been  substituted  into  eq  (5.19) 
and  the  resulting  expression  is  simplified,  it  becomes 

-  1  /s(2)   ^2[y^(t^^       At),   t^  -f  At] 

(S  (2)  \ 
—  +  ij  C2[(NR  -  DAy,   t^  +  At]  (5.24) 

2(2)  S  (2) 

-  1  +  S(2)   C2(NRAy,   t^  +  At)  =  ^  C2[.(NR  -  l)Ay,   t^  +  At]  , 

where 

2  D2At 

After  introducing  the  FORTRAN  notation  in  table  9  into  eq   (5.24),  it 
becomes 


1  +  S(2)    ^^^^^   TTMMPl)    +(^^Jjj-+  Ij   C2(NR,  TTMMPl) 

(5.26) 

S ( 2 )  ^  ( 7) 

-,    ,      '(^s    C2(NR  +  1,   TTMMPl)    =  C2  (NR,   TTMM)  . 


Table  9:     FORTRAN  Version  of  Quantities  in  Eq  (5.24). 

D2  =  DF2 

t       =  TTMM 
MM 

t.„,  +  At  =  TTMMPl 
MM 

At  =  DT 
Ay  =  DY 

2    •   DF2   •  03? 


R(2)  = 


DY** 


C2(Y    (t„^  +  At),   t       +  At)   =  CB(2) 
o     MM  MM 

C9((NR  -   l)Ay,    t       +  At)    =  C2 (NR,  TTMMPl) 
^  MM 

C2((NR  -  l)Ay,    t^^)    =  C2  (NR,  TTMM) 
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Equation  (5.26)  will  be  utilized  later  in  section  6.9  to  find  the  quan- 
tities PA(2),  QA(2)    such  that 


C2(NR,   TTMMPl)    =  PA(2)    CB(2)    +  QA(2)  (5.25') 

Schematic  summary  of  the  derivation  of  the  tridiagonal  FORTRAN  equation 
for  the  y-grid  point  NR, 

(5.19)  — >   (5.24)  — >  (5.26) 

Summary  note : 

Equation  (5.26)  is  used  in  section  6.9  to  determine  the  quantities 
PA(2),  QA(2)   associated  with  grid  point  NR. 

5.6.     The  Grid  Points  Lie  at  the  Moving  Boundary  and  NL  >  1    (The  Pre- 
liminary Equation) 

The  main  equation  associated  with  the  two  unknowns  Ci [z^iti^  +  At) ,  t^M 
+  At],  C2[yQ(tf4j4  +  At),  tf42^^  +  At]  located  at  the  grid  points  which  form 
the  interior  boundary  of  the  z-  and  y-regions,  respectively,  is  derived 
from  the  integral  form  of  the  conservation  of  number  boundary  condition 
eq   (2. 22)  ,  namely 

.z.(t,,,,+At)  ^y, 
^  '    ■     ■  —         -  C2(y,t)dy 


/o     MM  r^B 
Ci(z,t)dz  +  / 

^A  "^^o^W^t) 


(5.27) 


-  Di3  Ci 


+  Co 

y 


where  z     =    (NL  -  l)Az  and  y    =   (NR  -  l)Ay. 

Furthermore,   it  is  assumed  that  the  moving  boundary  z^{t)   is  such  that 
(NL  -  l)Az  ^  Zq  (t)   <_  NLAz,   for  all  times  t  between  and  tj^jyi  +  At. 

Similarly,   for  times  in  this  same  interval  the  moving  boundary  yQ(t)  is 
assumed  to  satisfy 

(NR  -  2)Ay  <_  y^(t)   <_  (NR  -  l)Ay  . 

With  these  restrictions,   eq   (5.27)   is  equivalent  to  the  differential 
boundary  condition  eq  (2.16). 

Before  proceeding  with  the  discretization  of  eq  (5.27),  it  is  convenient 
to  introduce  immediately  the  FORTRAN  notation  in  table  10. 
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Table  10:     FORTRAN  Symbols  Used  in  Eq  (5.28). 


Ci((NL  -   l)Az,   t...  +  At)   =  C1(NL,  TTMMPl) 

Ci({NL  -  l)Az,   t     )   =  CKNL,  TTMM) 
J-  MM 

C2((NR  -  l)Ay,   t..,.  +  At)   =  C2(NR,  TTMMPl) 
MM 

C2((NR  -  l)Ay,   t„)   =  C2  (NR,  TTMM) 
^  MM 

Ci  (z    (t       +  At)  ,   t       +  At)   =  CB(1,  TTMMPl) 
^     o     MM  MM 

Co(y    (t       +  At)  ,   t       +  At)   =  CB(2,  TTMMPl) 
o     MM  MM 

^l^^o^^MM^'   ^MM^   =  ™^ 


By  using  the  backward  difference  quotient  approximation  in  eq  (A7)  to 
approximate  the  time  derivatives  on  the  left-hand  side  of  eq  (5.27)  and 
using  the  trapezoidal  rule  to  approximate  the  integrals  in  the  result- 
ing expression,  we  find  (using  the  notation  introduced  in  table  10)  that 
the  left-hand  side  of  eq  (5.27)   is  equivalent  to 

S (1) Az  S (2) Ay 

  C1(NL,   TTMMPl)   +  CB(1,   TTMMPl)   +    CB(2,  TTMMPl) 

SM  (1) Az 

+  C2(NR,    TTMMPl)    -   Cl(NL,  TTMM) 

(5.28) 

+  CB(1,  TTMM)  -  ^^2At^'^  CB(2,  TTMM)  +  C2  (NR,  TTMM) 
+  higher  order  terms  in  Ay,  Az,  At  , 

where 

S(l)   -   iz   (t^,  +  At)   -    (NL  -  l)Az]/Az 
o  MM 

SM(1)   =    [z    (t     )   -    (NL  -  l)Az1/Az 
o  MM 

S(2)   -   [(NR  -  DAy  -  Y^it^  +  At)]/y  (5.29) 

SM(2)  =  [(NR  -  DAy  -  y   (t,,J]/Ay  . 

o  MM 

These  are  consistent  with  previous  definitions  of  S(l)   and  S(2), 

When  eqs    (5.14)   and   (5.22)    (re-expressed  in  the  FORTRAN  notation  of  ta- 
ble 10)   are  employed  to  approximate  the  spatial  derivatives  on  the 
right-hand  side  of  eq   (5.27)  and  these  approximations  along  with  eq 
(5.28)   are  substituted  into  eq   (5.27),   there  results  after  simplifica- 
tion and  up  to  terms  of  higher  order 
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^-^^C1(NL,   TTMMPl)    +  CB(1,    TTMMPl)    +  ^— — ^-^CB(2,  TTMMPl) 
R(l)  R(2) 

+  C2(NR,   TUVIMPl)    -  ^^,1]^    CKNL,   TTMM)    +  CB(1,    TTMM ) 

R  ( 1 ) 

-  ^   '  CB(2,   TTMM)    +  C2  (NR,  TTMM) 

R  ( 2 ) 

^  -Dl  I  -   I  Vs\^  )  ~   1'  TTMMPl) 

+    (S(l)    -   1)    C1(NL,   TTMMPl)  (5.30) 


where 


S(1)(A  s(i))  ™""'| 


*         -  S(2)(l  t  S(2))'^^''- 


 ^^^^    C2(NR,  TTMMPl) 

S  { 2 ) 

S  (2)  ) 

C2  (NR  +  1,    TTMMPDJ.  , 


1  +  S(2) 


RM 

^  =  m 

RM  =  ^  (5.31) 

Ay 

1 

DR  =  —  . 

D2 

In  order  to  eliminate  the  division  by  S(2)    (which  could  have  the  value 
zero)   in  the  last  term  on  the  right-hand  side  of  eq   (5.30),  we  eliminate 
that  term  by  substituting  from  the  relation 

1  1  -  S  (2) 

-  S(2)(l  +  S(2))   ^^^2'  ^       S(2)  ^2^^^' 

S  (2) 

+   -.    +  g(2)    C2(NR,  TTMMPl) 

(5.32) 

=  C2  (NR  +  1,    TTMMPl)    -  C2  (NR,  TTMMPl) 
1 


R(2) 


[-C2(NR,   TTMMPl)    +  C2  (NR,  TTMM)] 


35 


This  latter  equation  is  just  the  formula   [equivalent  before  simplifica- 
tion to  eq   (5.24)]   that  results  when  eqs   (5.20),    (5.21),   and   (5.22)  are 
substituted  into  eq   (5.19)   and  the  FORTRAN  notation  of  table  9  intro- 
duced . 


If  the  terms  in  the  brackets  on  the  right-hand  side  of  eq   (5.30)  are 
eliminated  by  substitution  from  eq   (5.32)   and  if  the  resulting  equation 
is  multiplied  through  by  S(l)    (to  prevent  possible  division  by  zero), 
then  we  obtain   [after  CB(2,  TTMMPl)   has  been  replaced  by   (m  =  SEG)  • 
CB(1,   TTMMPl)   using  boundary  condition  eq  (2.23)]: 


S(l)'- 
R(l) 


Dl 


SEG 


S(l)    •  S(2) 


1  +  S(l) 
S(l)2 


R(2) 

+    (1  -  S(l)) 


CB(1,  TTMMPl) 


R(l) 

R   •   S (1)    •   S (2) 

R(2) 


C1(NL,  TTMMPl) 


-  R  •  S(l) 


C2(NR,  TTMMPl) 


S(l)^    •  Dl 
1  +  S(l) 

S(l)    •  SM(1) 
R(l) 

R   •   SM(2)    •  SM(1) 


C1(NL  -   1,   TTMMPl)    +  R   •   S (1)    C2 (NR  +  1,  TTMMPl) 


R(2) 


C1(NL,    TTMM)    +  CB(1,  TTMM) 


CB(2,  TTMM) 


(5.33) 


+  S(l) 


R  •    SM(2)  R     \    ^  , 

+  — -T— |C2(NR,  TTMM) 


R(2) 


R(2)  / 


Further  FORTRAN  notation  in  table  11  can  be  introduced  into  eq  (5.33) 
and  by  using  intermediate  variables,   it  can  be  more  compactly  rewritten: 


DM1    •   CB(1,    TTMMPl)    =  DM2    •   C1(NL,   TTMMPl)    +  DM3   •    C2 (NR,  TTMMPl) 


(5.34) 


+  DM4 
where 


C1(NL  -   1,    TTMMPl)    +  DM5    •   C2 (NR  +   1,   TTMMPl)    +  DM6  , 


DM1  = 


DM2 


DM3  =^  - 


S(l) 


DFl 


SEG 


RT 


S(l)    •  S(2) 


R(l)  1  +  S(l) 

S(l)**2 


R(2) 


R(l) 

RT  •   S(l)    •  S(2) 


DM4  = 


DFl 


R(2) 
S  (1) **2 


+  DFl   •    (1  -  S(l))  , 
-  RT   •  S(l) 


1  + 


(2)  ) 


1  +  S(l) 


(5.34') 
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DM5  =  RT   •   S (1)    ,  and 


S(l)    •  SM(1) 

DM6  =   — r-   C1(NL,   TTMM)    +  CB  (1 ,  TTMM) 

R  ( 1 ; 

RT   •    SM(2)    •  S(l) 
+   TT^^  —  CB(2,  TTMM) 


+  S(l)    .  •  -  .   C2(NR,   TTMM)  , 

R  ( 2 ) 


and  where 

DF2   •  DZ 

R  =  RT  = 


R(l)  = 
R(2)  = 


DY 

2   •  DT 

. DZ**2  ' 

2    •   DF2    •  DT 


DY**2 

S(l)    =    (ZOPl  -    (NL  -   1)    •    DZ)/DZ   ,  (5.35) 

S(2)    =    ( (NR  -   1)DY  -  Y0P1)/DY  , 

SM(2)   =    (  (NR  -   DDY  -  y„(t     ))/DY   ,  and 

^  MM 

SM(1)   =    (z    (t     )    -    (NL  -  1)DZ)/DZ  . 
o  MM 

Table  11:     FORTRAN  Version  of  Symbols  in  Eq  (5.33). 


At  =  DT 

Ay  =  DY 

Az  =  DZ 

z    (t  + 
o  MM 

At) 

=  ZOPl 

y   (t  + 
-^o  MM 

At) 

-  YOPl 

Di   =  DFl 

D2  =  DF2 

SUMMARY 


Schematic  summary  of  the  derivation  of  the  FORTRAN  equation  for  the  z- 
grid  point  at  the  moving  boundary,   i.e.,  at   (^^(tj^  +  At),  when  NL  >  1; 

(5.27)  — >  (5.33)  — ^  (5.34)   and  (5.34') 
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Summary  notes : 


(1)  All  the  concentrations  on  the  right-hand  side  of  eq  (5.34), 
i.e.,   C1(NL,   TTMMPl) , . . . ,C2 (NR  +  1,   TTMMPl)   will  be  eliminated 
in  terms  of  known  quantities  and  thereby  CB(1,  TTMMPl)   will  be 
found.     This  is  carried  through  in  section  6.10. 

(2)  The  concentrations  in  DM6  are  known  because  they  are  evaluated 
at  time  TTMM. 

5.7.     The  Grid  Point  Is  at  the  Moving  Boundary  and  NL  =  1 

The  derivation  of  the  preliminary  equations  for  C 1  (  Zq  ( tj^j  +  At),   tf^  + 
At)    is  very  similar  to  the  derivation  of  the  preliminary  equation  for 
the  same  quantity  in  the  previous  section,   so  if  the  user  desires  more 
details  than  are  given  below,  he  should  refer  to  that  section. 

We  again  begin  with  the  conservation  eq   (5.27).     However,   instead  of 
approximating  the  spatial  derivative  at  z  =        by  eq  (5.14),  we  use 
boundary  condition  eq   (2.24).     Because  of  this  change,  we  find  in  place 
of  eq  (5.30) 

S  (1) 

— TT-T-   [C1(NL,    TTMMPl)    +  CB(1,   TTMI4P1)  ] 
R  ( 1 ) 

R  •  S(2) 


R(2) 

SM(1) 
R(l) 

R  •  SM(2) 


[CB(2,    TTMMPl)    +  C2(NR,  TTMMPl)] 


[C1(NL,    TTMM)    +   CB(1,  NR,  TTMM)] 

[CB(2,    TTMM)    +  C2(NR,   TTMM)]  (5.36) 


R(2) 


=  -    [+  C     Az(Cl(l,    TTMMPl)    -  C  )] 
p  out 


+  R 


1  1  -  S  (2) 

CB(2,   TTMi4Pl)    +   C2(NR,  TTMMPl) 


S(2)  (1  +  S(2)  )    ^    '    '  '  S(2) 

S(2) 


1  +  S(2) 


C2 (NR  +  1,  TTMMPl) 


When  the  term  in  the  second  parenthesis  on  the  right-hand  side  of  eq 
(5.35)   is  eliminated  by  using  eq   (5.32),  there  results 

S  (1) 

•     [C1(NL,    TTMMPl)    +  CB(1,  TTMMPl)] 
R  ( 1 ) 

+  ^  •   ^^^^    [CB(2,   TTMMPl)    +  C2(NR,  TTMMPl)] 
R(2) 
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^^777   [C1(NL,   TTMM)    +  CB(1,   TTMM) ] 
R  ( i ) 

R  •  SM(2) 


R(2) 


[CB(2,   TTMM)   +  C2(NR,   TTMM)]  (5.37) 


=  -  [C     Az(Cl(l,  TTMMPl)    -  C  J] 
p  out 


+  R    [C2(NR  +  1,   TTMMPl)    -  C2 (NR,  TTMMPl)] 
R 


R(2) 


[-C2(NR,   TTMMPl)    +  C2(NR,  TTMM)] 


By  rearranging  terms  and  introducing  the  FORTRAN  symbols  in  eq  (5.35) 
into  eq   (5.37),  we  get 

DM1    •   CB(1,    TTMMPl)    =  DM2    •   Cl (NL,   TTMMPl)    +  DM3    •   C2 (NR,  TTMMPl) 

+  DM4  •  C1(NL  -  1,  TTMMPl)  (5.38) 
+  DM5    •   C2(NR  +  1,   TTMMPl)    +  DM6  , 

where 

S(l)        RT   '3(2)    •  SEG 

DMl  =    +    , 

R(l)  R(2) 

DM2  =  -  ^-TTT"  +  CONPRO   •   DZ  , 
R(l) 

=  -  iif  -  ^(^ '  ^) '  ^'-'^'^ 

DM4  =  0  , 

DM5  =  RT   ,  and 

DM6  =  ^^11^     [C1(NL,   TTMM)    +  CB(1,  TTMM)] 
R  ( X ) 

■  •    SM(21fc3(2,  TTMM)] 


R(2) 

RT   •   SM(2)    -  RT 


[C2 (NR,   TTMM) ] 


R(2) 

+  DZ   •   CONPRO   •   CONOUT  , 

and  where  the  FORTRAN  symbols  above  are  the  same  as  in  eq  (5.35)  with 
the  addition  of  the  ones  in  table  12. 
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Table  12:     Additional  FORTRAN  Symbols  for  Eq  (5.38). 


C     =  CONPRO 
P 

C     ^  -  CONOUT 
out 


Notice  that  the  quantities  DM3,  DM5,   and  DM6   [with  DZ  •  CONPRO  •  CONOUT 
added  to  it  in  eq   (5.34')]   are  equal  to  their  counterparts  in  eq  (5.38') 
when  S (1)  -  1.     Hence,  by  setting  S(l)   to  1  when  we  are  computing  these 
quantities,   i.e.,   in  the  case  NL  =  1,  we  can  use  the  same  FORTRAN  expres 
sion  as  is  used  to  compute  these  quantities  in  eq   (5.34')   when  NL  >  1. 

Of  course  in  computing  in  this  way  we  must  leave  the  S(l)   in  DMl  and 
DM2  of  eq   (5.38')   unchanged.     This  can  be  done,  when  NL  =  1,  by  intro- 
ducing the  quantity  SAVE  which  is  equal  to  S(l)   and  replacing  S(l)  where 
it  appears  in  DMl  and  DM2  by  the  symbol  SAVE. 

SUMMARY 

Schematic  summary  for  obtaining  the  preliminary  equation  for  the  z-grid 
point  at  the  moving  boundary  when  NL  1. 

(5.27)  — >   (5.37)  --^  (5.38)   and   (5.38')  . 

Summary  notes : 

(1)  Since  the  coefficients  of  the  concentrations  on  the  right- 
hand  side  of  eq   (5.38)   are  the  same  as  those  on  the  right- 
hand  side  of  eq   (5-34),  the  concentrations  from  both  of  these 
preliminary  equations  can  be  eliminated  by  the  same  procedure 
which  will  be  described  in  section  G.IO. 

(2)  The  actual  computation  of  the  quantities  in  eq   (5.38')  are 
carried  through  as  described  above  in  computing  block  #4  of 
subroutine  CENTER,  see  section  7.2. 

5.8.     The  Index  of  the  Z-Grid  Point  is  1  and  NL  >  1 

The  derivation  of  the  tridiagonal  algebraic  equation  that  is  associated 
with  the  z-grid  point  1   (when  the  moving  boundary  has  reached  the  sec- 
ond z-grid  point)  begins  with  the  conservation  of  number  eq  (2.20) 

(5.39) 

where  z^^  =  0  and  Zg  =  Az/2  denote  the  positions  of  the  flux  boundaries. 


d 

dt 


/ 


Ci  (z, t)dz  = 


-Di3  Ci 
z  ^ 


+  Di9  Ci 


z=z. 


A 
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After  approximating  the  time  derivative  on  the  left-hand  side  of  eq 
(5.38)  by  the  back  difference  quotient  approximation  eq   (A7)   and  the 
integrals  by  quadrature  formulas  with  a  node  at  the  left-hand  end  point, 
we  find 


d 

dt 


2   (t  +At) 
o  MM 


Ci(z,t)dz  =  — —  [Ci(0  •  Az,   t      +  At)   -  Ci  (0   •   Az,   t  )] 

^  2At       ^  MM  ^  MM 


Az 


(5.40) 


+  0(Az^  +  At^/At)  . 
By  substituting  from  eq  (2.24), 


C   (C1(0  •   Az,  t.„,  +  At)   -  C  ) 
p  MM  out 


for  Di 9  Ci 
z 


on  the  right-hand  side  of  eq  (5.39)  and  then 


z=z. 


D,  (Ci  (Az,   t      +  At)   -  Ci  (0  •  Az,   t      +  At) ) 
1     1  MM  ^  MM 


for  DiS  Ci 
z 


on  the  right-hand  side  of  eq  (5.39),  that  equation 


z=z 


B 


yields  together  with  eq  (5.40)  after  simplification 


Az^ 


+  Az   •   C     +  Di 
P 


Ci(0  •  Az,   t,,,  +  At)   =  DiCi(Az,   t^,^  +  At) 
^  MM  Mil 


+  :   +  Az   •   C     •   C  ^ 

2  •  At  p  out 

Az2 


(5.41) 


After  introducing  into  eq  (5.41)  the  FORTRAN  notation  in  table  13,  that 
equation  becomes 


/2dt\  ^ 
\DZ2j 


+  DZ   •   CONPRO  +  DFl 
Did,  TTMM) 


CKy,   TTMMPl)    =  DFl    •   C1(2,  TTMMPl) 


(5.42) 


2DT 
DZ**2 


+  DZ   •   CONPRO  •  CONOUT 
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Table  13:     FORTRAN  Version  of  Quantities  in  Eq  (5.41). 


At  =  DT 

Az  =  DZ 

C     =  CONPRO 
P 

C         =  CONOUT 
out 

Di   =  DFl 

Ci(0   •   Az,   t  „  +  At)   =  Cl(l,  TTMMPl) 
MM 

Ci(0    •    Az,    t,.)    -  Cl(l,  TTMM) 
J-  MM 

Co(Az,   t       +  At)   =  C2{2,  TTMMPl) 

^  MM 


Equation   (5.42)   is  the  FORTRAN  tridiagonal  equation  which  will  give  rise 
to  the  "P,  Q"  quantities  associated  with  the  z-grid  point  1  and  unknown 
Cl(l,  TTMMPl).     This  is  carried  through  in  section  6.4. 

SUMMARY 

Schematic  summary  of  the  derivation  of  the  FORTRAN  tridiagonal  equation 
for  the  z-grid  point  1  when  the  moving  boundary  is  beyond  z-grid  point  2 , 

(5.39)  — >   (5.41)  — >  (5.42). 

Summary  note : 

The  FORTRAN  eq   (5.42)   leads  to  values  of  Pl(l),   Ql(l)   that  are  com- 
puted in  computational  block  #24  of  the  main  program.     See  section 
5.4  for  the  derivation  of  these  "P,  Q"  quantities. 

5.9.     The  Index  of  the  Z-Grid  Point  is  1  and  NL  =  1 

This  equation  is  different  from  all  the  others  in  that  it  is  not  devel- 
oped from  an  integral  conservation  equation.     It  is  basically  derived 
from  combining  eq   (2.14)   with  a  Taylor  series  expansion. 

From  a  spatial  Taylor's  formula  expansion  of  Ci (z,   tMM  +  At),  we  have 


Ci  (z   (t      +  At)  )   -  Ci  (0,   t      +  At)   +  9  Ci 
^     o    MM  ^  MM  z  ^ 


z   (t      +  At) 
o  MM 

z=0  (5.43) 


+  1/2  aci 


(z    (t^^  +  At))2  +  0(Az3)  . 
o  MM 


z=0 
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By  defining 


S(l)   =  z   (t      +  At)/Az 
o  MM 


and  substituting  this  together  with   (from  eq  (2.14)) 

into  eq  (5,43),  it  becomes  after  rearranging  and  up  to  higher  order  terms 

(S(l)Az) 


DnCi (z    (t       +  At) ,   t       +  At)   =  DnCi (0,   t       +  At)   +  Di 9  Ci 
i^oMM  MM  ^^  MM  -^z^ 


z=0 


+  1/2  3^Ci 


(S(l)    •  Az)^ 


(5.44) 


z=0 


When  the  backward  difference  approximation, 

3^Ci  =   (Ci(0,   t^^  +  At)   -  Ci(0,   t,^))/At  +  0(At) 
t  .       MM  MM 

and  the  expression  from  eq  (2.24), 


Dia^ci 


=  C    (Ci  (0,   t^^  +  At)   -  C     ,  )  , 
p  MM  out 


(5.45) 


z=0 


are  substituted  into  eq  (5.44),  it  becomes,  after  rearrangement  (ignoring 
higher  order  terms) , 


At  •   2  p 


Ci(0,  t      +  At) 
MM 


=  DiCi(z    (t^  +  At),   t^^  +  At) 
o     MM  MM 


(5.46) 


+  C,  (0,   t  ) 
^  MM 


(S(l)    ♦  Az) 
2At 


+  S (1)    •   Az   •   C     •  C 

p  out 


By  transferring  the  symbols  in  eq  (5.46)  through  the  correspondences  in 
table  14,   that  equation  becomes 

(S(l)    •  DZ)**2 

DFl  +  ^       2   '   DT    ^^"^^    *  *  '  TTMMPl) 


=  DFl  CB(1,  TTMMPl)  +  Cl(l,  TTMM) 
+  S(l)    •   DZ   •   CONPRO   •   CONOUT  . 


(S(l)    •  DZ)**2 
2DT 


(5.47) 
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Table  14:     FORTRAN  Version  of  Quantities  in  Eq  (5.46). 


Di  =  DFl 
Az  =  DZ 
At  =  DT 

Cl(0,  +  At)   =  Cl(l,  TTMMPl) 

MM 

Cl(0,    tj    =  Cl(l,  TTMM) 
MM 

Ci (z    (t       +  At) ,   t       +  At)   =  CB(1,  TTMMPl) 
^     o     MM  MM 

C     =  CONPRO 
P 

C     ^  =  CONOUT 
out 


FORTRAN  tridiagonal  eq   (5.47)   will  lead  to  the  quantities  Pl(l),  QKD 
in  section  6.3.     These  are  the  "P,  Q"  associated  with  z-grid  point  1  when 
the  moving  boundary  has  not  yet  reached  the  second  z-grid  point. 

Notice  that  when  t^j^^  =  0  the  quantity  Cl(l,0)   is  unknown,   i.e.,   it  is  not 
given  as  data  in  the  problem.     Thus,   the  equation  we  use  in  place  of  eq 
(5.45)   for  the  first  iteration  MM  =  1  when  tj,^  =  0  is 

(DFl  +  S(l)    •   DZ   •   CONPRO)    •   Cl(l,   TTMMPl)    =  DFl   •   CB(1,  TTMMPl) 

(5.48) 

+  S(l)    •   DZ   •   CONPRO   •    CONOUT  . 

This  is  the  equation  we  would  have  arrived  at  if  we  had  truncated  the 
Taylor's  formula  in  eq   (5.43)   to  first  order  terms  and  proceeded  just  as 
above,  and  eq  (5.48)   is  consistent  with,  i.e.,  the  same  as,   the  formula 
we  would  get  by  approximating  eq  (5.45)  directly. 

SUMMARY 

Schematic  summary  of  the  derivation  of  the  FORTRAN  tridiagonal  equation 
associated  with  the  z-grid  point  1  when  NL  =  1, 

(5.43)  — >  (5.45)  — ^  (5.47)   or   (5.48)    (use   (5.48)   when  t^  =  0). 

Summary  note : 

The  "P,  Q"  terms  that  result  from  eq   (5.47)   or  eq   (5.48)   are  de- 
rived in  sections  6.2  and  6.3. 
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6.     METHOD  OF  SOLVING  THE  ALGEBRAIC  EQUATIONS 


6.1.     The  Derivation  of  the  "Double  Sweep"  "P,  Q"  Coefficients  Associated 
with  the  Grid  Points 

In  this  section  we  put  each  of  the  tridiagonal  FORTRAN  equations  associ- 
ated with  the  fixed  grid  points  in  "double  sweep"  form.     Specifically,  we' 
manipulate  these  tridiagonal  FORTRAN  eqs   (3.2)   and   (3.3)   into  the  forms 

C1(J,  TTMMPl)  -  P1(J)  CKJ  +  1,  TTMMPl)  +  Ql  ( J)  (6.1) 
for  J  =  1,...,NL  -  1,  and 

CKNL,  TTMMPl)  =  PA(1)  CB(1,  TTMMPl)  +  QA(1)  .  (6.2) 
Also,   at  the  fixed  grid  points  in  the  y-region, 

C2(J,  TTMMPl)  =  P2(J)  CI (J  +  1,  TTMMPl)  +  Q2 ( J)  (6.3) 
for  J  =  NR  +  1,...,N2,  and 

C2(NR,   TTMMPl)    =  PA(2)    CB(2,   TTMMPl)    +  QA(2)    .  (6.4) 

Once  again  the  specific  composition  of  the  double  sweep  "P,  Q"  factors 
associated  with  a  grid  point  depends  on  the  location  of  the  grid  point 
and  the  position  of  the  moving  boundary  at  time  TTMMPl.     The  following 
table  15  completely  classifies  all  the  cases  that  arise.     Of  course,  be- 
cause the  "P,  Q"  coefficients  for  a  grid  point  are  obtained  from  the  tri- 
diagonal equation  associated  with  that  grid  point,  this  table  is  the 
same  as  table  3  except  for  Case  I (e)  and  Case  II  (g)  which  were  not  in- 
cluded in  table  3   (since  they  are  so  trivial) . 


Table  15:     Guide  to  the  Subsections  Where  the  "P,  Q"  Coefficients  As- 
sociated With  Each  of  the  Fixed  Grid  Points  Are  Derived  and 
Described. 


Case  I :         The  index  NL  equals  1 : 

(a)  The  first  z-grid  point,   see  sections  6.2  and  6.3. 

(b)  The  index  of  the  y-grid  point  is  greater  than  NR,   see  sec- 
tion 6.8. 

(c)  The  index  of  the  y-grid  point  is  NR,  see  section  5.9. 

(d)  The  grid  point  lies  at  the  moving  boundary,   see  section 
6. 10. 

(e)  The  largest  y-grid  point,  which  has  index  N2  and  is  lo- 
cated at  y  =  BDRY  =  NYO  •   DY,   see  section  6.7. 
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Case  II:       The  index  NL  is  greater  than  1. 


(a)  The  index  of  the  z-grid  point  is  1,  see  section  6.4. 

(b)  The  index  of  the  z-grid  points  ranges  from  2  to  NL  -  1, 
see  section  6.8. 

(c)  The  index  of  the  z-grid  point  is  NL,   see  section  6.6. 

(d)  The  grid  point  is  at  the  moving  boundary,  see  section  6.10. 

(e)  The  index  of  the  y-grid  point  is  greater  than  NR,  see  sec- 
tion 6.8. 

(f)  The  index  of  the  y-grid  point  is  NR,  see  section  6.9. 

(g)  The  index  of  the  y-grid  point  is  N2  and  is  at  y  =  BDRY, 
see  section  6.7. 


6.2.  The  First  Z-Grid  Point  When  NL  =  1  and  t      =  TTMM  =  0 

MM 

After  dividing  both  sides  of  eq   (5.48)  by   (DFl  +  S(l)    •  DZ  •  CONRPO) , 
we  get 

Cl(l,   TTMMPl)   =  Pl(l)    •   CB(1,   TTMMPl)   +  Ql(l)  (6.5) 

where 

Pl(l)   =  DM  •  DFl 

Ql(l)   =  DM  •    (Cl(l)    •   DM3  +  DM1   •   CONPRO  •  CONOUT) 
DM1  =  S(l)    •  DZ 
DM2  =  0 
DM3  =  DM2 

DM  =    (DFl  +  DM2  +  DMl   •   CONPRO)"^  . 

The  introduction  of  the  zero  quantities  DM2,  DM3  into  the  above  equa- 
tions enables  us  to  use  below  in  section  6.3  the  same  expression  DM  but 
with  DM2,   DM3  ^  0. 

SUMMARY 

Equations   (6.5')  are  calculated  in  computational  block  #23  of  the  main 
program   (see  sec.  7.1). 

6.3.  The  First  Z-Grid  Point  When  NL  =  1  and  t      =  TTMM  >  0 

MM 

f S (1)    •   DZ) **2 

When  eq  (5.47)   is  divided  through  by  DFl  +   ^   ^     2  DT    ^  ^^^^    *  ' 

CONPRO,  there  results 

Cl(l,    TTMMPl)    =  Pl(l)    CB(1,   TTMMPl)    +  Ql (1) 

with 


(6.5') 


(6.6) 
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(6.6') 


Pl(l)    -  DM   •  DFl 

QUI)    =  DM   •    (Cl(l,   TTMM)    •   DM3  +  DMl   •   CONPRO   •  CONOUT) 
DM  =    (DFl  +  DM2  +  DMl   •  CONPRO) 
DMl  =  S(l)    •  DZ 
DM2  =   (DMl) 2/2   •  DT 
DM3  =  DM2  . 

SUMMARY 

Note  that  all  quantities  in  eq  (6.6')  are  the  same  as  in  eq  (6,5')  except 
for  DM2  and  DM3.     The  calculation  of  the  quantities  in  eq   (6.6')  takes 
place  in  computational  block  #23  of  the  main  program  (see  sec.  7.1). 

6.4.  The  Index  of  the  Z-Grid  Point  is  1  and  NL  >  1 

(DZ)  **2 

After  dividing  eq   (5.42)   by    +  DZ  •  CONPRO  +  DFl),   it  becomes 

^     ^  2   •  DT 

Cl(l,   TTMMPl)   =  Pl(l)    •   Cl(2,   TTMMPl)   +  Ql  (1)  (6.7) 

where 

Pl(l)    =  DFl   •  DM 

DZ**2 

01(1)    =  DM  Cl(l,    TTMM)    •   +  DZ   •   CONPRO  •   CONOUT  (6.7') 

2   •  DT 

DM  =  ^  £_  +  DZ   •    CONPRO  +  DFl  . 

2    •  DT 

SUMMARY 

The  quantities  in  eq   (5.7')   are  computed  in  computational  block  #24  of 
the  main  program   (see  sec.  7.1). 

6.5,  The  Index  of  the  Z-Grid  Point  is  Between  1  and  NL  -  1  and  NL  >  2 

We  begin  with  eq  (5.7)  which  holds  for  the  indices  J  -  2,...,NL  -  1.  As- 
sume that  for  one  of  the  indices  in  this  range 

C1(J  -   1,   TTMMPl)    =  P1(J  -   1)    •    CKJ,   TTMMPl)    +  Ql  ( J  -   1)    .  (6.8) 

After  substituting  this  value  of  Cl ( J  -  1,  TTMMPl)   into  the  left-hand 
I  side  of  eq  (5.7)  and  rearranging  the  terms,  there  results 

CKJ,   TTMMPl)    =  PI  (J)    •   C1(J  +  1,   TTMMPl)    +  Ql  (J)  (6.9) 

where 
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PI 
Ql 

Dl 


(J) 
(J) 
(J) 


=  -A2R/D1(J) 

=  (AF(J)  -  QKJ  -  1)A2L)/D1(J) 
=  A2M  =  PI (J  -  1)    •   A2L  . 


(6.9') 


Notice  that  eq  (6.8)  with  J  -  2  holds  because  of  eq  (6.7).  Hence,  by  in- 
duction we  get  eqs    (6.9)   and   (6.9')   for  J  =  2,...,NL  -  1. 


Summary  notes : 

(1)  The  arrays  Pi (J) ,  Dl  (J)  are  computed  in  computational  block 
#3  of  subroutine  TRIDNL. 

(2)  The  quantities  A2L,...,A1M   (see  table  5)   are  computed  in  sub- 
routine INTER  which  is  called  by  subroutine  TRIDNL  when  Pi  (J)  , 
Dl(J)  are  being  computed. 

(3)  The  arrays  AF(J),  Ql  (J)   are  calculated  in  computational  block 
#25  of  the  main  program. 

6.6     The  Index  of  the  Z-Grid  Point  is  NL  and  NL  >  1 

When  eq  (5.17)   is  multiplied  through  by  S(l)  (since  S(l)  may  be  zero, 

this  prevents  possible  division  by  zero)  and  the  terms  are  transposed, 
there  results 


SUMMARY 


DM1  C1(NL,   TTMMPl)    =  DM2  CB(1,  TTMMPl) 


+  DM3  C1(NL  -   1,   TTMMPl)    +  DM4 


(6.10) 


with 


DM  = 


1  +  S(l) 


DM  = 


DM2 


DF1/(1  +  S(l))    =  DFl/DM 


(6.10' ) 


DM4 


DM3 


DFl  •  S  (1)/(1  +  S  (1)  )  =  DFl 
S  (1) 

;.,  :    Cl    (NL,    TTMM)  . 


DM2 


We  also  have  from  eq  (6.9)  with  J  =  NL  -  1, 


C1(NL  -  1,    TTMMPl)    =  P(l)    C1(NL,    TTMMPl)    +  Q(l) 


(6.11) 


where  we  have  used  the  abbreviations 


P(l)   =  P1(NL  -  1) 


Q(l)   =  QKNL  -  1) 


6.11') 
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When  eq  (6.11)  is  substituted  into  the  right-hand  side  of  eq  (6.10)  and 
the  result  is  simplified,  we  find 


C1(NL,   TTMMPl)    =  PA(1)    CB(1,   TTMMPl)    +  QA(1)  (6.12) 

where 

PA(1)    =  DM2/ (DM1  -  DM3    •  P(l)) 

(6.12') 

QA(1)   =    (DM4  +  Q(l)    •   DM3) /(DM1  -  DM3   •  P(l))  . 

SUMMARY 

The  quantities  in  eqs  (5.10')  and  (6.12')  are  evaluated  in  computational 
I   block  #3  of  subroutine  CENTER. 

:i    6.7.     The  Index  of  the  Y-Grid  Point  is  N2  =  NYO  +  1  and  is  located  at 
y  =  BDRY 

I  The  boundary  condition  for  C2 (   ,   )   at  the  boundary  point  y  =  BDRY  is 
j    (see  eq   (2. 12) ) 

C2(N2,   TTMMPl)   =  CBULK  .  (6.13) 

j    We  rewrite  this  simple  equation  for  use  later  on  in  the  more  convenient 
l'  form 

C2(N2,   TTMMPl)   ^  P2  (N2)   C2 (N2  -  1,   TTMMPl)   =  Q2(N2)    ,  (6.14) 
if  we  assume   (and  we  do) 

P2(N2)   =  0 
Q2(N2)   =  CBULK  . 

SUMMARY 

||  The  quantities  in  eq   (6.14^)  are  calculated  in  computational  block  #13 
of  the  main  program.     These  quantities  are  used  in  block  #26  of  the  main 
:.rogram  (see  sec.   7.1)   and  block  #4  of  siabroutine  INTER  (see  sec.  7.3). 


6.8.     The  Index  of  the  Y-Grid  Point  is  Between  NR  and  N2 


We  start  with  tridiagonal  eq  (5.10)  which  holds  for  J  =  NR  +  1,..., 
112  -  1. 


Assume  for  one  of  the  indices  in  this  range  that 

C2(J  +  1,   TTMMPl)   =  P2(J  +  1)   C2(J,  TTMMPl)   +  Q2(J  +  1)    .  (6.15) 

After  substituting  this  value  of  C2 (J  +  1,  TTMMPl)   into  the  left-hand 
side  of  eq  (5.10),  with  J  in  place  of  J  -  1,  and  rearranging  terms,  it 
becomes  for  J  =  NR  +  1, . . . ,  N2  -  1 
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C2{J,  TTMMPl) 


=  P2(J)   C2(J  -  1,   TTMMPl)    +  Q2(J) 


(6.16) 


where 

1 


P2(J)   =  - 


D2  (J) 

Q2(J)  =  (BF{J)  -  Q2(J  +  1)  B2R)/D2(J)  (6.16') 
D2(J)   =  B2M  +  P2(J  +  1)B2R 

or,  since  B2R  =  1  from  table  7, 

D2 (J)   =  B2M  +  P2 (J  +  1)  . 

Notice  that  eq  (6.15)  with  J  =  N2  -  1  holds  because  of  eq  (6.14).  Hence, 
by  induction  we  get  eq   (6.16)  and  eq  (6.16')   for  J  =  NR  +  1,...,N2  -  1. 

SUMMARY 

Summary  notes : 

(1)  The  arrays  P2 (J) ,  D2(J)  are  calculated  in  computational  block 
#4  of  subroutine  TRIDNL. 

(2)  The  quantities  B2M,  BlM   (the  other  B's  equal  1)   are  computed  in 
subroutine  INTER  which  is  called  by  subroutine  TRIDNL  when 

P2  (J) ,  D2 (J)   are  computed. 

(3)  The  arrays  Q2(J),  BF(J)  are  calculated  in  computational  block 
#26  of  the  main  program. 

6.9.     The  Y-Grid  Point  has  Index  NR 

After  transposing  terms  in  the  tridiagonal  FORTRAN  eq  (5.26),  it  takes 
the  form 

DM1  C2(NR,   TTMMPl)    =  DM2  CB(2,  TTMMPl) 

(6.17) 

+  DM3  C2(NR  +  1,   TTMMPl)   +  DM4 


where 


DMl  =  1  + 


DM2 


R(2) 
1 


^  +  '^'2>  (6.17') 
DM3  =  S  (2)/(l  +  S  (2)  ) 

S  (2) 

DM4  =  ^l^j    C2 (NR,   TTMM)  . 


We  also  have  from  eq  (6.15)  with  J  =  NR  +  1 
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C2  (NR  +  1,   TTMMPl)    =  P(2)    •   C2  (NR,   TTMMPl)    +  Q(2) 


(6.18) 


where  we  have  used  the  abbreviations 

P(2)=P2WR+1) 
Q(2)    =  Q2(NR  +  1)  . 

When  eq   (5.18)   is  substituted  into  the  right-hand  side  of  eq   (6.17)  and 
the  result  is  simplified,  we  find 

C2(NR,   TTMMPl)   =  PA(2)   CB(2,   TTMMPl)   +  QA(2)  (6.19) 

where 


PA (2)   =  DM2/ (DM1  -  DM3   •  P(2)) 

QA(2)   =    (DM4  +  Q(2)    •   DM3) /(DMl  -  DM3   •  P(2)) 

SUMMARY 


(6.19') 


The  quantities  eqs   (6.1?')   and   (6.19')   are  evaluated  in  computational 
block  #3  of  subroutine  CENTER. 

6.10.     The  Determination  of  CB(1,  TTMMPl)   in  Terms  of  Known  Quantities 

The  values  of  the  "P,  Q"  coefficients  in  the  eqs   (6.1)   through  (6.4)  are 
determined  by  the  coefficients  in  the  tridiagonal  eqs   (3.2)   and  (3.3) 
and  these  in  turn  are  composed  of  known  mesh  parameters,  material  parame- 
ters and  concentrations  evaluated  at  time  TTMM.     Hence,  by  using  eqs  (6.1) 
through   (6.4),  we  can  eliminate  the  concentrations  evaluated  at  time 
TTMMPl  that  occur  on  the  right-hand  side  of  the  "preliminary"  eq  (5.34) 
and  thereby  determine  the  value  of  CB(1,  TTMMPl)  entirely  in  known  quan- 
tities . 

More  precisely,  substituting  eqs  (6.11)  and  (6.18)  into  the  right-hand 
side  of  the  "preliminary  equation,"  i.e.,  eq  (5.34)  for  NL  >  1  (and  eq 
(5.38)   for  NL  =  1),  we  get 


(6.20) 


DMl  CB(1,   TTMMPl)    =  DM8   •   Cl (NL,  TTMMPl) 

+  DM9  C2(NR,  TTMMPl)   +  DMlO 

where 

DM8  =  DM2  +  DM4    •  P(l) 

DM9  =  DM3  +  DM5   •   P(2)  (6. 20' ) 

DMlO  =  DM6  +  Q(l)    •   DM4  +  Q(2)    •   DM5  . 

Finally  eliminating  the  concentrations  at  the  grid  point  NL  and  NR  on  the 
right-hand  side  of  eq   (6.20)   by  use  of  eqs    (6.12)   and   (6.19)   leads  to  the 
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expression  for  CB(1,  TTMMPl)   in  terms  of  known  quantities 

DMll 

CB{1,   TTMMPl)    =  ,  (6.21) 

with 

DM12  =  DM1  -  DM8   •   PA(1)    -  SEG   •   DM9   •   PA (2) 
DMll  =  DM10  +  QA(1)    •   DM8  +  QA{2)    •   DM9  , 


and  where 


DM1 
DM2 

DM3 

DM4 
DM5 
DM6 


.     eq   (5.39)    if  NL  =  1 

are  given  m         /c  o/in    •     ,.tt  i 
eq   (5.  34)   if  NL  >  1 


and  furthermore. 


PA(2) 
QA(2) 


are  given  in  eq  (6.19') 


^^2)  given  in  eq  (6.18  ) 

PA(1)  .          .  ,„/. 

are  given  m  eq  (6.12  ) 

QA(1) 


P(l) 
Q(l) 


are  given  in  eq  (6.1l') 


(6.21') 


(6.22) 


R(l)     is  given  in  eq  (5.16  ) 
R(2)     is  given  in  eq  (5.35) 
RT  =  R  is  given  in  eq  (5.31) 

SUMMARY 

The  quantities  in  eqs  (6.21),  (6.2l')  and  all  the  quantities  listed  below 
eq  (5.2l')  are  computed  in  computational  block  #4  of  subroutine  CENTER  to 
find  the  value  of  CB(1,  TTMMPl). 
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7.      DETAILED  DESCRIPTION  OF  THE  PROGRAM 


7,1.     Block-by-Block  Detailed  Description  of  the  Calculational  Procedure 
in  the  Main  Program 

As  is  shown  in  the  program  listing  (Chapter  11) ,  the  overall  calculations 
in  DISTRIB:     VERSION  1  are  broken  down  into  a  set  of  consecutive  blocks. 
The  computations  in  each  of  these  blocks  have  a  unified  purpose  which  is 
described  with  pertinent  comments  for  each  of  these  blocks  #1  through 
#36  in  this  chapter.     In  Chapter  8,  the  logic  in  these  blocks  is  flow 
charted.     The  precise  calculations  carried  through  in  each  of  these 
blocks  are  shown  in  the  program  listing  in  Chapter  11,  and  the  material 
in  the  following  blocks  is  intended  to  illuminate  the  contents  by  orga- 
nizing and  elaborating  on  the  material  in  the  listing  rather  than  to  re- 
describe  the  contents  in  a  different  language. 

IMPORTANT :     Quantities  with  names  which  appear  in  the  program  listing  v 
but  do  not  appear  in  the  glossary  are  not  needed  or  used  in  DISTRIB: 
VERSION  1  and  therefore  should  be  ignored,  i.e.,  all  material  in  the 
listing  that  is  not  in  a  block  or  glossary. 

Block  #1:     Input  Declaration,   Type  Statements,  Dimensionalize  Arrays 
and  Common  Blocks 

The  arrays  SRL{     )   and  YRL(     )   are  declared  real  because  of  certain  con- 
straints in  the  plotting  routines  where  these  arrays  are  employed.  The 
fact  that  they  are  dimensionalized  by  200  means  that  only  200  data  points 
can  be  plotted  in  a  single  call  of  the  plot  routine. 

Arrays  that  have  an  element  associated  with  each  grid  point,   e.g.,  C2  ( J) 
(see  glossary)   are  dimensionalized  by  2000.     This  means  the  maximum  num- 
ber of  grid  points  is  2000  and  the  maximum  number  of  mesh  widths  in  ei- 
ther region  must  be  less  than  2000. 

Since  the  number  of  cells  of  the  array  Cl (     )  actually  used  is,  at  time 
TTMMPl,   equal  to  NR  •  NZO,   the  input  quantities  NZO  and  TFINAL  (which 
together  with  the  equation  of  the  moving  boundary  determine  the  largest 
possible  value  NR  takes)   must  be  constrained  so  that  the  product  NR  • 
NZO  is  always  less  than  the  dimension  of  Cl (  ). 

Many  of  the  arrays  are  dimensionalized  in  common  blocks  because  they  are 
also  used  in  subroutines. 

Block  #2:     Input  the  Proper  Problem  and  Material  Parameters  with  Accom- 
panying Format  Declarations 

The  definition  of  each  of  the  input  quantities  is  to  be  found  in  the 
glossary,  where  references  are  given  to  places  in  the  main  body  of  this 
documentation  which  contain  more  elaborate  descriptions  of  these  quan- 
tities.    The  list  of  input  and  output  is  given  in  section  9.2  in  a  sys- 
tematic fashion  and  in  classified  form. 
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Real  input  is  read  under  format  4  and  real  output  printed  under  format 
5.     Integer  input  and  output  are  under  format  1. 

A  table  listing  the  input  data  is  found  in  section  9.2  along  with  ranges 
of  values  for  the  input  parameters  that  have  worked  well  in  trial  runs. 
The  units  for  the  input  data  are  micrometers,  hours.     The  problem  with 
CONOUT  =  0  is  homogeneous  in  the  concentrations,  and  therefore  the  in- 
put concentrations  are  normalized  to  CBULK  =1.     By  homogeneous  is  meant 

the  value  of  C. (x,t)   for  CBULK  =  A  is  equal  to  C. (x,t)    [for  CBULK  =  12 
1  1 

txmes  A. 

Block  #3:     Store  Input  in  Terms  of  Their  Original  Dimension  Units 

In  the  next  block  certain  input  parameters  are  dimensionalized  (rescaled) 
Consequently,   these  input  quantities  are  stored  in  their  original  units 
in  the  array  PARAM  so  they  can  be  printed  out  later  in  block  #34  in  their 
original  units.     Also,  before  its  defining  quantities  are  changed,  the 
quantity  SCALE2    (in  micrometers  which  is  used  to  scale  the  distaTice  along 
the  y-coordinate  system,   i.e.,   in  the  silicon,  in  the  plots  that  are 
printed  in  block  #34  is  calculated  and  stored  for  later  output.     The  rea- 
son SCALE2  does  not  appear  in  the  plot  routine  is  that  its  role  there  is 
played  by  UN2  which  represents  the  same  physical  quantity  as  SCALE2  but 
is  in  units  of  BLTH  rather  than  micrometers.     The  fact  that  DY,  YOPl, 
etc.,  when  they  appear  in  the  plot  routines  are  in  BLTH  units  requires 
the  use  of  UN2  instead  of  SCALE2.     A  detailed  description  of  the  dimen- 
sionalization  procedure  is  given  in  section  4.3. 

Block  #4:     Redimension  the  Input  Data  by  Introducing  Computationally 
Efficient  Units 

The  dimensionalization  carried  out  here  is  described  in  section  4.3  in 
detail.     Note  that,   in  Version  1,  DFC  is  set  equal  to  1  in  a  declare 
card.     Thus,  during  the  calculations,   DF2  will  have  the  value  1  and  DFl 
will  be  the  ratio  DF1/DF2.     The  computational  formula  for  BTM  is  given 
in  eq   (4. 10) . 

Block  #5:     Compute  the  Values  of  the  Grid  Mesh  Widths 

This  is  the  step  described  in  section  4.1  where  DY  and  DZ,   the  mesh 
widths  in  the  y-  and  z-regions,   respectively,  are  calculated.     The  com- 
putational formula  for  DY  is  eq   (4.1)   and  eq   (4.3)   for  DZ. 

Block  #6:     Compute  the  Value  of  ZCRIT  =  CNSTX. 

In  this  step  CNSTX,  which  is  another  name  for  ZCRIT,   is  calculated,  see 
section  4.2.     The  quantities  CNA  and  CNB  are  just  two  intermediate  vari- 
ables used  in  computing  DNSTX  and  ZOPl,   Z0P2  later  on. 
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Block  #7:     Compute  the  Time  Increment  DT  for  the  First  Iteration  of  the 
Main  Loop  78 

In  this  block  the  main  purpose  is  to  compute  the  initial  value  of  At 
(=  DT) .     Recall,  from  section  4.2,  that  in  the  "wet  oxidation"  Case  I 
the  time  TLIN  =  CNTAU  =  x  =  0,   so  the  motion  of  the  moving  boundary  is 
always  parabolic.     In  this  case  DT  is  defined  as  the  time  it  takes  the 
moving  boundary  to  move  the  distance  DM  =  WFRC  •  DZ,   having  started  at 
Z  =  0.     The  number  WFRC  must  be  positive,  real,  and  less  than  1,  for 
the  program  to  be  stable.     In  the  case  of  a  "dry  oxidation,"  Case  II  in 
section  4.2  which  starts  at  Z  =  0    (in  VERSION  1)   TLIN  ='CNTAU  =  t  >  0, 
the  boundary  moves  linearly  for  the  time  period  CNTAU  with  a  speed 
SLOPE  =  ZCRIT/CNTAU.     Hence,   the  time  DT  the  moving  boundary  needs  to 
transverse  the  distance  DM   (where  DM  is  as  above)   is  DM/SLOPE.  Values 
of  WFRC  that  have  been  successfully  tested  range  between  0.01  and  0.33. 

We  have  required  TLIN  to  be  less  than  or  equal  to  10~^  rather  than  zero 
in  order  to  distinguish  between  the  wet  and  dry  oxidation  cases.  This 
does  not  seem  now  to  be  necessary,  but  we  have  left  this  as  it  is  be- 
cause in  all  practical  cases  TLIN  is  always  much  greater  than  10~^  when- 
ever it  is  greater  than  zero,   so  no  harm  is  done. 

The  quantity  DTS  is  used  to  store  the  quantity  DT  for  recall  at  a  later 
step,  see  computational  block  #17. 

When  CNTAU  =  t  =  0,   the  quantities  ZCRIT  =  CNSTX  and  SLOPE  which  are  not 
used  but  still  appear  in  later  formulas  are  assigned  the  values  0  and  1, 
respectively,  which  do  not  foul  up  these  formulas,  for  example  by  caus- 
ing division  by  zero. 

In  case  x  >  0   (i.e. ,  dry  oxidation) ,  WFRC  must  certainly  be  chosen  so 
that  DM  <  ZCRIT  for  the  above  prescription  for  DT  to  make  sense.  In 
all  the  cases  we  have  run,  this  has  been  naturally  satisfied  when  data 
in  the  literature  are  satisfied. 

Block  #8:     Computation  and  Storage  of  Initial  Value  DT  in  Hours 

The  quantity  DT  just  calculated  in  the  previous  block  is  multiplied  by 
BTM   (the  unit  in  which  time  has  been  rescaled)   to  determine  its  magni- 
tude in  units  of  hours.     It  is  stored  in  the  array  PARAM  for  printout 
in  block  #34. 

Block  #9:     Calculation  of  Initial  Values  of  NL,  NR  and  Value  of  N2  for 
the  First  Iteration  of  the  Main  Loop  78 

The  quantity  Nl  is  another  name  for  the  index,  NL,  of  the  first  grid 
point  with  a  z-coordinate  less  than  or  equal  to  the  z-coordinate  of  the 
moving  boundary  at  the  beginning  time  of  the  main  iteration  loop  78. 
Because  YOl ,  ZOl  are  zero  in  Version  1,  Nl  has  the  value  1,  as  it  should 
since  at  this  point  in  the  program  the  moving  boundary  is  at  z  =  0  (re- 
call we  are  computing  here  NL  for  the  first  iteration  of  loop  78)  NR 
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(which  is  the  index  of  the. first  grid  point  that  has  a  y-coordinate 
greater  than  the  y-coordinate  of  the  moving  boundary  at  the  start  of 
the  iterations  in  loop  78)  has  the  value  2, 

The  quantity  N2  is  the  index  of  the  last  grid  point  on  the  y-axis  and 
since  NYO  is  the  number  of  grid  mesh  widths  in  the  interval   (0,  BDRY) 
(see  below  eq  (2.11)   for  the  definition  of  BDRY)  N2  =  NYO  +  1   (see  com- 
putational formula  eq  (4.2)). 

Block  #10:     Calculation  of  Initial  Values  of  C2  (  ) 

The  initial  values  of  array  C2 (     )  at  all  the  grid  points  are  determined 
from  the  prescription  in  eq  (2.13)   of  C2(y,0).     Because  of  this  prescrip- 
tion, the  grid  points  from  NR  through  NR  +  NH   (NH  is  an  input  integer) 
are  assigned  the  value  CMAX   (input  value)  while  all  the  following  grid 
points  are  assigned  the  value  CBULK  (input  value) . 

Since  in  VERSION  1  no  oxide  is  initially  present,  there  is  no  necessity 
to  assign  array  Cl  (     ) . 

Notice  that  C2 (1)   is  essentially  assigned  through  the  assignment  of 
CB(2).     This  is  because  CB(2)  denotes  by  definition  the  concentration 
in  the  silicon  at  the  position  of  the  moving  boundary  and  initially  the 
moving  boundary  is  located  at  the  position  of  the  first  y-grid  point. 

Block  #11:     Calculation  of  the  Value  of  TTMM  for  MM  =  1 

This  is  the  step  in  which  TTMM,  the  time  at  the  beginning  of  the  MMth 
iteration  of  loop  78,  for  the  first  iteration  is  calculated.     Since  ZOl 
in  VERSION  1  has  been  assigned  the  value  zero,  in  a-  declaration  card, 
this  quantity  obviously  has  the  correct  value  which  is  zero. 

Block  #12:     Calculation  of  the  Initial  Values  of  SM(  ) 

The  quantities  SM(1),  SM(2)   are,  respectively,  the  distance  between  grid 
point  NL  and  ZOPl   (the  position  of  the  moving  boundary  in  the  z- 
coordinate  system  at  time  TTMMPl)   and  the  grid  point  NR  (for  more  de- 
tails see  figs.   5  and  6) . 

The  quantity  FNL  represents  the  niamber  of  mesh  widths  between  z-grid 
point  1  and  z-grid  point  NL,  and  FNR  is  the  number  of  mesh  widths  be- 
tween y-grid  point  NR  and  y-grid  point  N2   (the  last  grid  point  on  the 
y-axis  in  the  domain  of  the  silicon) . 

Since  ZOl  and  YOl  are  declared  zero  in  VERSION  1,  SM(1)  =  0  and  SM(2)  = 
1  as  they  should  be  when  we  are  calculating  these  quantities  for  the 
initial  iteration  of  loop  78. 
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Block  #13:     Calculation  of  the  "P,  Q"  Coefficients  for  Grid  Point  N2 

As  has  been  explained  in  section  6.7,  the  satisfaction  of  the  boundary 
condition  C2(N2)   =  CBULK   (see  eq   (2.12))   at  y  =  BDRY  is  equivalent  (in 
the  way  this  problem  is  being  solved)   to  assuming  P2(N2),  Q2(N2)  are 
assigned  in  such  a  way  that  C2(N2)   =  P2(N2)   C2  (N2  -  1)   +  Q2(N2)  [see 
eqs   (6.14)   and   (6.14') ] . 

The  values  of  P2(N2),  Q2(N2)   are  used  in  subroutine  TRIDNL. 

Block  #14:     Beginning  of  the  Main  Loop  which  Calculates  Concentrations 
at  Time  =  TTMMPl 

At  the  beginning  of  this  loop,  78,  the  quantities  in  the  arrays  CI  (  ), 
C2 (     )   representing  the  concentrations  in  the  z-  and  y-regions  at  each 
grid  point  have  their  values  at  time  TTMM.     At  the  end  of  the  MMth  itera 
tion,  these  same  quantities  will  have  their  proper  values  at  time  TTMMPl 

Before  the  main  loop  begins,   certain  starting  values  for  quantities  used 
in  the  loop  are  set.     The  definitions  and  use  of  these  quantities  ZOPS, 
TJ,  etc.,  are  explained  in  later  steps. 

The  overall  manner  in  which  the  new  concentrations  are  calculated  is 
explained  in  section  3. 

Block  #15:     Setting  the  Print-out  Control  Parameter 

The  integer  KW  controls  print-out  in  block  #30.     This  data  will  be 
printed  out  after  the  first  iteration,  then  skip  KPRINT  -  1  iterations 
of  loop  and  print  out  the  data  again,  and  so  on. 

Block  #16:     Increasing  the  Distance  that  the  Moving  Boundary  Travels 
in  an  Iteration  of  Loop  78 

During  the  MMth  iteration  of  loop  78,  the  quantity  ZOACEL  measures  that 
fraction  of  the  mesh  width  DZ  that  the  moving  boundary  transversed  dur- 
ing the  previous  iteration  of  loop  78.      (If  MM  -  1,  ZOACEL  =  0  since 
ZOPl  =  0.   and  ZOPS  =  0.     In  this  case,   see  block  #14.)     If  ZOACEL  is 
less  than  0.2  and  KPOW   (defined  below)   is  greater  than  100,  then  the 
magnitude  to  DT  is  increased  by  0.05  in  order  to  accelerate  the  distance 
the  boundary  moves  during  the  MMth  iteration  of  loop  78.     This  accelera- 
tion is  put  into  the  algorithm  in  order  to  decrease  the  number  of  itera- 
tions that  are  required  in  order  to  compute  the  solution  to  a  problem  in 
which  the  moving  boundary  moves  a  certain  specified  distance. 

Of  course,   subroutine  TRIDNL,  which  computes  and  stores  quantities  that 
depend  on  DT,  has  to  be  recomputed  and  the  new  value  of  DT  is  placed  in 
DTS.     When  leaving  the  loop  in  which  DT  has  been  increased,   the  value  of 
the  control  parameter  KPOW  is  increased  by  100  so  that  DT  will  not  be 
increased  again  for  at  least  another  100  iterations  of  loop  78. 
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The  quantities  0.2  and  100  mentioned  above  are,  of  course,  somewhat 
arbitrary,  but  the  above  values  have  worked  satisfactorily  in  a  great 
number  of  trial  runs. 

Block  #17:     Adjustments  Required  When  Moving  Boundary  Reached  New  Z-Grid 
Point  on  Previous  Iteration  of  Loop  78. 

The  quantity  TJ  will  be  computed  in  a  later  step,  block  #20,  in  such  a 
way  that  when  and  only  when  the  moving  boundary  reached  a  new  z-grid 
point  during  the  previous  step  will  TJ  have  a  value  greater  than  or 
equal  to  1.     Hence,  assuming  T J  >_  1 ,  NL,  which  is  the  first  grid  point 
less  than  or  equal  to  the  position  of  the  moving  boundary  at  the  time 
when  MMth  iteration  of  loop  78  starts,  has  to  be  increased  in  value  by 
one.     The  quantity  FNL   (the  number  of  mesh  widths  to  the  left  of  NL) 
also  must  be  increased  by  1  and  NLMl  =  NL  -  1,  NLM2  =  NL  -  2,  increased 
accordingly. 

Also,  CI  at  the  new  grid  point,  NL,  must  be  defined. 

Furthermore,   since  DT  has  been  adjusted   (see  block  #20)  during  the  pre- 
vious iteration  of  loop  78,  this  quantity  is  reset  to  the  value  it  had 
during  the   (MM  -  2)th  iteration,   specifically  DTS.     Also,  TRIDNL,  the 
subroutine  computing  the  arrays  PI,  Ql,  etc.,  depends  on  DT   (besides  a 
new  PI,  Dl  being  added  in  the  y-region)   and  hence  needs  to  be  recomputed. 

The  value  SM(1)    (see  sees.   5.6  and  5.7  for  its  meaning)  must  be  set 
equal  to  zero  because  during  the   (MM  -  l)th  iteration  the  moving  boun- 
dary reached  a  new  grid  point  in  the  z-region. 

All  these  changes  are  effected  when  TJ  >  1.     Of  course,   in  the  contrary 
case  the  implication  is  that  the  moving  boundary  had  not  yet  reached  a 
new  z-grid  point  and  the  above  changes  are  not  required. 

Block  #18:     Computing  the  Time  and  Position  of  the  Moving  Boundaries  at 
the  End  of  the  MMth  Iteration  of  Loop  78 

In  this  step  we  are  simply  computing  what  the  value  of  the  time  will  be 
at  the  end  of  the  MMth  iteration  of  loop  78.     Since  at  the  beginning  of 
the  loop,  t  has  the  value  TTMM  by  definition  and  since  time  is  made  to 
increase  by  DT  during  the  loop,  the  time  at  the  end  of  the  loop  TTMMPl 
is  given  by  TTMM  +  DT.     After  this  value  of  TTMMPl  is  computed,  the  lo- 
cation of  the  moving  boundary  in  the  z-  and  y-regions,  respectively, 
ZOPl,  YOPl  can  be  computed  through  the  use  of  formulas  in  section  4.2 
and  eq   (2.8).     Remember,   ZW  and  TW  equal  zero  in  VERSION  1. 

Block  #19-     Adjustments  Required  When  Moving  Boundary  Reached  New  Y- 
Grid  Point  on  Previous  Iteration  of  Loop  78 

The  quantity  SES  is  the  fraction  which  the  signed  distance  between  NR 
and  YOPl  is  of  the  mesh  width  DY  =  AY.     If  YOPl  is  beyond  the  y-grid 
point  NR  (during  the  MMth  iteration  of  loop  78) ,  then  SES  is  negative 
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and  NR,  FNR,   and  NRPl  =  NR  +  1  have  to  be  increased  by  1.     When  SES  is 
negative  at  the  beginning  of  the  MMth  loop,  it  means   (see  the  next  block) 
the  moving  boundary  is  located  at  time  TTMM  on  a  grid  point,   and  hence 
SM(2)   should  be  set  equal  to  1. 

Block  #20:     If  the  Moving  Boundary  Moves  for  a  Time  2DT,  Will  It  Reach 
a  New  Z-Grid  Point? 

During  the  MMth  iteration  of  loop  78,  we  look  at  where       (   )   is  at  time 
TTMMP2  =  TTMMPl  +  DT.     This  position  is  denoted  as  Z0P2.     Then  we  com- 
pute TJ  which  is  that  fraction  the  signed  distance  between  grid  point 
NL  and  ZOPl  is  of  the  mesh  width  AZ  =  DZ .     If  this  fraction  is  less 
than  one,  it  means  that  during  the   (MM  +  l)th  iteration  of  loop  78  the 
moving  boundary  will  not  reach  a  new  z-grid  point,   and  we  proceed  with 
the  MMth  iteration.     On  the  other  hand,   if  TJ  is  greater  than  or  equal 
to  1,  we  know  that  the  moving  boundary  will  reach  or  cross  grid  point 
NL  +  1  on  the  next  iteration   (MI'I  +  1)  .     Therefore,  we  decide  to  change 
the  value  of  DT  so  that  while  still  in  the  MMth  iteration  of  loop  78, 
Zq  travels  all  the  way  to  the  next  grid  point.     At  the  time  TTr4MPl  when 
this  happens,    i.e.,   Zo(TTr4MPl)   =   (FNL  +1),  DZ  is  calculated   (as  shown 
in  the  listing)    for  the  two  cases  CNTAU  =  0  or  CNTAU  >  0.  (Remember 
that  ZW  =  TW  =  0  in  VERSION  1.)     Finally,  we  calculate  the  new  value  of 
DT  that  will  take  the  moving  boundary  Zq  (   ),   starting  at  time  TTMM,  to 
its  new  position.    .Lastly,  we  recompute  TRIDNL  because  the  arrays  Pi, 
P2,  etc.,   computed  there  depend  on  the  value  of  DT. 

Block  #21:     Calculation  of  the  Parameters  S(l),   S(2)   Needed  for  Subrou- 
tine CENTER 

The  quantities  FNL,  FNR,   etc.,  have  now  been  reset   (if  it  was  necessary). 
These  quantities  S.(l)   and  S(2)    (see  sees.   5.4  to  5.7  for  details)  are 
used  in  computing  the  coefficients  of  the  unknowns  in  the  algebraic 
equations  that  will  be  solved  later  in  subroutine  CENTER. 

Block  #22:     For  the  First  Iteration  of  Loop  78,  Call  TRIDNL 

In  the  first  iteration  of  loop  78,   the  arrays  P2,  D2  that  are  needed 
later  and  which  are  calculated  in  subroutine  TRIDNL  have  not  been  cal- 
culated before  in  the  program. 

Notice  that  whenever  a  new  NL  occurs,  and  hence  the  quantities  in  TRIDNL 
need  to  be  updated,   this  is  done  in  step  #20. 

Block  #23:     Calculating  of  Pl(l),  QKD  When  the  Moving  Boundary  Zq 
Has  not  Yet  Reached  Grid  Point  Two 

In  the  case  NL  =  1,  the  moving  boundary  at  time  TTMMPl  has  still  not 
reached  the,  second  z-grid  point.     Therefore,  Pl(l),  Ql(l)  must  be  cal- 
culated using  formulas  in  sections  6.2  and  5.3,  specifically  eqs  (6.5') 
when  TTim.  =  0  and   (6.6')   when  TTMM  >  0. 


Of  course,  DM,  DMl,   etc.,   are  just  intermediate  variables. 
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For  the  computational  significance  of  Pl(l),  Ql  (1)  ,   see  again  sections 
5.2  and  6.3.     Pl(l)   is  not  calculated  in  subroutine  TRIDNL  like  the  rest 
of  its  elements,  and  Ql(l)   is  calculated  here  rather  than  with  the  other 
elements  of  Ql  in  block  #25. 

Block  #24:     Computation  of  Pl(l),  Ql(l)  When  Zq  Has  Moved  Beyond  the 
First  Grid  Mesh  Widths 

In  this  case   (NL  >  1) ,   the  moving  boundary  has  traveled  through  the 
first  z-mesh  width  and  Pl(l),  Ql(l)   are  to  be  calculated  by  the  formulas 
in  section  6.4,   specifically  eq   (6.7').     Pl(l)   is  not  calculated  in  sub- 
routine TRIDNL  like  the  rest  of  its  elements,   and  Ql(l)   is  calculated 
here  rather  than  with  the  other  elements  of  Ql  in  block  #25. 

Block  #25:     Calculation  of  Arrays  AF  (  ),  Ql  (   )   VJhen  There  are  Regular 
Points  in  the  Z~Region 

When  NL  >  2,   there  are  grid  points  with  cell  widths   (see  sec.   5.1)  of 
magnitude  DZ  in  the  z-region,  and  the  arrays  AF  and  Ql  associated  with 
these  points  are  therefore  computed.      [This  amounts  to  saying  that  the 
arrays  AF (  ),  Ql (  )   that  are  associated  with  these  points  are  calculated 
by  using  the  formulas  in  table  5.2  and  eq   (6.9').  ] 

E*or  the  significance  of  the  array  Ql  in  the  computational  procedure, 
see  chapter  3.     The  array  AF  plays  the  role  of  the  quantity  Zi+  in  eq 
(3.2);   see  also  eq   (5.7).     Array  AF  is  used  to  compute  array  Dl  in  sub- 
routine TRIDNL.     The  array  Ql,  which  is  composed  of  AF  and  Dl,   is  used 
in  block  #29. 

Block  #26:     Calculation  of  Arrays  BF  and  Q2  Associated  with  the  Grid 
Points  in  the  Y-Region 

The  array  BF  is  used  to  compute  Q2  later  in  this  block. 

For  further  elucidation  of  the  computational  significance  of  these  ar- 
rays,  see  the  references  to  BF  in  section  5.3,   table  7,  and  eq  (5.10). 
See  section  6.8  and  eqs   (6.16)  and   (6.16')   for  the  significance  and 
definition  of  Q2(J).     Notice  P2 (N2) ,  Q2(N2)   are  calculated  in  block 
#13.     Arrays  P2,  Q2  are  used  in  block  #28. 

Block  #27:     Calculation  of  CB(1),  CB(2),  CI (NL) ,  C2 (NR)   at  Time  TTMMPl 

The  quantities  above  at  the  end  of  the  MMth  iteration  of  loop  78  (i.e., 
at  time  TTMMPl)   are  computed  in  subroutine  CENTER  which  is  documented 
in  detailed  block  form  in  section  7.2. 

Block  #28:     Computation  of  C2 (   )   at  Time  TTMMPl 

By  making  use  of  C2 (NR)   that  was  computed  in  subroutine  CENTER  in  block 
#27,  we  can  carry  out  the  backward  sweep;   see  steps  #2  and  #3  in  chapter 
3  that  generate  C2(J)   for  J  =  NR  +  1  to  N2  -  1. 
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Block  #29:     Computation  of  Array  CI (   )   at  Time  TTMMPl 

In  this  case  NL  =  1,  the  value  of  Cl(l),  which  is  the  only  grid  point 
besides  CB(1)   in  the  y-region,  is  already  known,  since  it  has  been  com- 
puted in  subroutine  CENTER  as  CI  (NL) .     However,  when  NL  is  greater  than 
one,   the  calculation  of  CI  (J)   for  J  =  1,  NL  -  1  must  be  completed  by  a 
backward  sweep,  as  described  in  section  3.1,  steps  #2  and  #3. 

Block  #30:     Print-out  Concentrations  When  KW  =  0 

The  integer  KW  is  zero  when   (MM  -  1)   is  divisible  by  KPRINT, 

The  concentrations  are  printed  under  a  time  TTMMPl  which  is  in  units  of 
BTM  (see  sec.   4.3).     The  user  is  advised  not  to  print  out  every  concen- 
tration in  each  region  because  there  are  so  many,  especially  in  the  y- 
region  (i.e.,  the  silicon). 

Also  printed  out  are  the  values  of  the  concentrations  at  the  boundary 
and  CI (NL) .      (Note  that  C2 (NR)  will  always  be  printed  with  the  other  con- 
centrations in  the  y-region. ) 

Block  #31:     Reset  Quantities  in  Preparation  for  Next  Iteration  of  Loop 
78 

Unless  a  new  grid  point  is  reached  by  the  moving  boundary  in  the  next 
iteration,  the  new  values  of  SM(1),  SM(2),  TTMM  are  as  given   (see  figs. 
5  and  6) .     If  in  the  next  iteration  of  loop  78  a  new  grid  point  is 
reached,   the  values  are  corrected  appropriately  at  that  time,  in  blocks 
#17  and  #19. 

Block  #32:     Test  to  Determine  Whether  TTMMPl  Has  Reached  TFINAL 

The  quantity  TFINAL   (INPUT)   is  the  time  we  want  to  stop  the  problem 
and  TIMEND  is  the  value  of  TFINAL  in  units  of  BTM  which  is  the  unit  of 
TTMMPl.     Hence,   as  long  as  TIMEND  is  greater  than  TTMMPl,  we  proceed  to 
the  next  iteration  of  loop  78.     If  MME  is  not  large  enough  to  reach 
TFINAL,  a  diagnostic  to  this  effect  is  printed  before  control  shifts  to 
statement  79  and  then  the  end  card. 

Block  #33:     Shift  Program  to  Statement  800  in  VERSION  1 

The  time  is  greater  than  or  equal  to  TFINAL  and  in  VERSION  1  the  value 
of  JSTOP  is  set  equal  to  1  in  a  data  card.     Hence,  control  sends  the 
program  to  statement  800  where  printing  and  plotting  occur  prior  to  the 
end  of  the  program. 

Block  #34:     Print-out  of  Data  and  Plot  of  Concentrations  at  the  Final 
Value  of  TTMMPl 

The  detailed  description  of  the  plotting  subroutine  is  given  in  Appendix 
2  and  section  4.4. 
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Block  #35:     Redimensionalize  Parameters  to  Original  Units 

Parameters  such  as  ZOPl  and  TTMMPl   (at  the  last  iteration  of  loop  78 
that  occurs)  are  changed  into  micrometer  and  hour  units  and  printed  out 
with  other  pertinent  parameters. 

Block  #36:     Control  in  VERSION  1  Sends  Program  to  Its  End 

The  integers  JSTOP  and  JSAVE  have  the  values  one  and  zero   (see  below 
data  cards)   in  VERSION  1;   consequently,   in  this  case  the  program  is  sent 
to  the  end. 

If  TIMEND  has  not  been  reached  by  TTI-lMPl  before  the  number  of  itera- 
tions MME  of  loop  78  is  completed,  a  diagnostic  message  to  this  effect 
is  printed  out. 

Block  #37:     This  Block  Contains  Statement  79  and  Ends  Program. 
7.2.     Block-by-Block  Description  of  Subroutine  CENTER 

This  subroutine,  which  is  called  in  computational  block  #27  of  the  main 
program,   is  where  the  quantities  CB(1),  CI (NL) ,  CB(2),  C2 (NR)   are  evalu- 
ated at  time  TTMMPl. 

Block  #0.     Insert  Input  Through  COMMON 

The  following  input  for  the  subroutine  passes  in  through  COMMON:     S(  ), 
SM(    ),   N2M1,   N2,   NLMl,   NL,   NR,   CI (    ),   C2 (    ),    CB (    ),   DFl,   DF2 ,  CONPRO, 
CONOUT,    DT,   DZ,   DY,   PI,   P2,   Ql,   Q2 . 

The  output  is  fed  to  the  main  program  through  COMMON. 
Block  #1:     Computation  of  the  Quantities  P(2),  Q(2) 

These  quantities  defined  in  eqs   (5.18)   and   (6.18')   are  used  to  deter- 
mine PA(2),  QA(2)   in  computational  block  #3. 

Block  #2:     Compute  P(l),  Q(l) 

These  quantities   [see  eqs    (6.11)   and   (6.11^)]   are  used  when  NL  >  1  to 
determine  PA(1),   QA(1).     When  NL       1,   the  role  played  by  PA(1),  QA(1) 
in  linking  the  value  of  CB(1)   to  the  value  of  Cl(l)   is  played  by  Pl(l), 
Ql(l),  which  are  computed  in  block  #23  of  the  main  program.  Therefore, 
at  the  end  of  block  #3,  we  set  PA(1)   =  Pl(l),  QA(1)   =  QKD   when  NL  -  1. 

Block  #3:     Computation  of  PA(L),  QA(L)    for  L  =  1,2 

The  computation  formula  for  these  quantities  when  L  =  1  is  described 
in  eqs    (6.10'),    (6.1l'),  and   (6.12').      [See  eqs   (5.16')   and   (5.25)  for 
the  definition  of  R(l)   and  R(2);  also  see  tables  7  and  8.]     The  compu- 
tation formulas  of  the  above  quantities  in  the  case  L  =  2  are  described 
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in  eqs  (6.17'),  (6.18'),  and  (6.19').  The  significance  of  PA(L),  QA(L) 
is  seen  from  eqs   (5.26')  and   (5.18)   or  eqs   (6.12)   and  (6.19). 

Block  #4:     Computation  of  CB(1),  CB(2)  at  Time  TTMMPl  in  Terms  of  Knovm 
Quantities 

The  concentration  at  the  moving  boundary  in  the  z-region  CB(1,  TTMMPl) 
CB(1)   is  computed  by  using  formula  eqs   (6.21)   and   (6.2l')  and  the  aux- 
iliary eq  (6.22).     To  find  CB(2)  =  CB(2,  TTMMPl),  we  just  use  the  boun- 
dary condition  in  eq   C2.23).     For  the  significance  and  explanation  of 
how  SAVE  is  used,   see  the  end  of  section  5,7,  before  SUMMARY. 

Block  #5:     Calculation  of  CI (NL) ,  C2 (NR)   at  Time  TTMMPl 

We  introduce  the  abbreviations 

C(l)  =  C1(NL,  TTMMPl)  =  CI  (NL) 
C(2)   =  C2(NR,   TTMMPl)   =  C2  (NR) 

and  compute  these  desired  quantities  from  CB(1),  CB(2),  PA(L),  QA(L) 
known  from  blocks  #1  to  #4  by  employing  eqs   (6.12)   and  (6.19). 

7.3.     Block-by-Block  Description  of  Subroutine  TRIDNL 

The  input  for  this  subroutine  is  the  output  of  subroutine  INTER,  DT, 
DZ,  DY,  DFl,  DF2,  and  CONPRO.  The  arrays  Pi,  Dl,  P2,  D2  in  contrast 
to  the  Ql,  Q2  arrays  do  not  depend  on  concentrations.  The  output  of 
this  subroutine  is  used  both  in  the  main  program  blocks  #25,  #26,  #28, 
and  #29  and  subroutine  CENTER  (see  blocks  #1  and  #2) .  The  output  of 
the  subroutine  is  PI (J),  Dl(J)  for  J  =  2,  NL  -  1  and  P2  (J) ,  D2(J)  for 
J  =  NR  +  1,    . . . ,  N2  -  1. 

Block  #1:     Call  Subroutine  INTER 

By  calling  this  subroutine,  we  input  A2L,  AIM,  B2L,  BIM 

through  COI-IMON.     These  quantities  are  used  in  blocks  #3  and  #4  here. 

Block  #2:     Calculation  of  Pl(l)   when  NL  >  1 

When  NL  >  2,  this  value  of  Pl(l)   is  used  in  the  next  block.  Calcula- 
tion formulas  of  the  above  quantities  are  found  in  eqs   (6.7)  and  (6.7') 
(This  is  a  duplication  of  Pl(l)   calculated  in  block  #24  of  the  main  pro 
gram  and  is  not  necessary  in  VERSION  1.) 

Block  #3:     Calculation  of  Arrays  Pi (J) ,  Ql(J) 

The  calculation  formulas  for  these  quantities  when  J  =  2,  NL  -  1  is 

eq   (6.9').     These  quantities  are  used  in  blocks  #25  and  #29  of  the  main 

program. 
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Block  #4:     Calculation  of  Arrays  P2  (J) ,  D2(J) 

The  calculation  formulas  for  these  arrays  when  J  =  NR  +  1,   . . . ,  N2  -  1 
are  given  in  eq  (6.16').     Notice  that  P2 (N2)   is  calculated  in  block  #13 
of  the  main  program  and  used  here. 

7.4.     Description  of  Subroutine  INTER 

This  subroutine  is  so  trivial  we  only  describe  it  briefly. 

The  routine  uses  the  value  of  DT,  DZ,  DY,  DFl,  DF2  to  compute  the  value 
of  A2L,   A2M,  A2R,   B2L,   B2M,   B2R,  AIM,   BlM  for  use  in  subroutine  TRIDNL 
(see  blocks  #3  and  #4) . 

Since  the  values  of  A2L,  BlM  depend  on  DT,  every  time  this  quan- 

tity changes  these  quantities  need  to  be  recomputed.     This  explains  why 
subroutine  TRIDNL  is  called  whenever  DT  changes  in  the  main  program, 
i.e.,  because  TRIDNL  calls  INTER  (see  block  #1  in  TRIDNL)   and  thereby 
corrects  the  values  of  A2L,  BlM  for  the  new  DT  and,  consequently, 

simultaneously  corrects  the  values  of  the  arrays  PI,  Dl,  P2,  D2  for  the 
values  of  the  new  DT. 
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8.     FLOW  DIAGRAMS 
8.1.     Flow  Diagram  of  Main  Program  in  DISTRIB 


#1 


#2 


#3 


#4 


#5 


^        START  ^ 


r 

DECLARATIONS , 
DIMENSION,  AND 
COM>'ON  BLOCK- 
STATEMENTS 

r 

INPUT  DATA 
AND  FORI^iAT 
STATEMENTS 

f 

STORE  INPUT 
IN  ORIGINAL 
PHYSICAL  UNITS 

f 

REDIMENSION 
THE  INPUT 
DATA 

COMPUTE  THE 
MAGNITUDES  OF 
THE  GRID  MESH 
V;iDTHS 

A 


65 


0 


#6 


#7 


#8 


COMPUTE  THE 
VALUE  OF  CNSTX 
=  ZCRIT 


COMPUTE  DT  FOR 
THE  FIRST 
ITERATION  OF  THE 
MAIN  LOOP  78 


COIJVERT  DT 
TO  HOUR  UNIT 
AND  STORE 


CALCULATION  OF  THE 
VALUES  OF  GBZD 
INDICES  NL,  NR, 
AND.  N2  FOR  THE 
FIRST  ITERATION  OF 
THE  MAIN  LOOP  78 


© 


#10 


f 

SET 
NH  + 

NH  = 
NR 

f 

DO  96, 

J  ---  KR,  NH 
C2  (J)    =  CMAX 

DO  97, 
J  =  NH 
C2(J)  = 

+  1,  N2 

CBULK 

r 

SET  CB(2)  = 
CMAX 

1 

f 

67 


CALCUTJ^TE  TTMM 
FOR  FIRST  ITER- 
ATION OF  THE 
MAIN  LOOP  78 


#12 


CALCULATION  OF  SM(1) 
SM(2)    FOR  THE  FIRST 
ITERATION  OF  THE 
MAIN  LOOP  78 


#13 


#14 


CALCULATION  OF 
THE  "P,Q"  CO- 
EFFICIENTS FOR 
GRID  POINT  N2 


BEGIN  MAIN  LOOP 
78  TO  FIND  CON- 
CENTRATIONS AT 
TIME  TTMMPl, 
SET  CONSTANTS  IN 
LOOP 


i 


68 


#15 


#16 


SET  PRINa-OUT 
PAFii^METER  ia-J 


T 


COMPUTE  DT 
TIME  /vCCELER- 
ATION  SIGNAL 
ZOACEL 


YES 


NO 

r 

UPVALUE 

SIGNAL 

KPOV;,  KPOW  = 

KPCW  + 

1 

<^  KPOW 

>  100 

NO 


YES 


INCREASE  DT. 
CALL  SUBROU- 
TINE TRIDNL. 
DL-VALUE  KPOW  BY 
100. 


RESET  20PS  -4 


69 


#17 


▼ 


UPVALUE  NL  AKD; 
FNL  BY  1.  SET 
CI  (ML)    =  CB(1) 


#18 


CAiyjUI^TE 

NLMl       NI,  - 

1 

4  

NLH2  =  NL  - 

2 

NO 


YES 


SET  DT  =  DTt 
CALL  TRIDNL 
SET  SM(1)  = 


CALCUIJ^TE  TTMIIPI 
FOR  tmth  ITERATION 
OF  xMAIN  LOOP  78  AND 
.POSITION  OF  nOVING 
BOUNDARIES  AT  TTMMPl 


V 

0 

70 


0 


#19 


CALCULATE 
THE  VALUE  OF 
SES  FOR  t-lMth 
ITERATION  OF 
LOOP  78 


YES 


SET  NEV?  VALUES 
OF  NR,  JNR  AND 
SMi2) 


SET  VALUES  OF 
NR  +  1,  NR  + 
2  AND  NRMl 


#20 


CALCUIJ^TE 
QUA-NTITIES 
NEEDED  FOR 
EVALUATING 
TJ 


#21 


1 

CALCUIiATx-- 
VALUE  OF 
TJ 

YES 


•CALCULATE  NEW 
VALUES  OF  ZOPl, 
YOPl,   DT  AND 
CALL  TRIDNL 


CALCULATE  VAL- 
UES OF  S (1)  , 
S  (2)   FOR  tmth 
ITERATION  OF  , 
MAIN  LOOP  78  ' 


72 


© 


#22 


CALL  TRIDNL 
FOR  USE  IN 
THE  FIRST 
ITEPJ\TION  OF 
LOOP  78 


#23 


#24 


<r  NL 

>   1  ^ 

NO 

f 

CALCULATE 

P]  (NL)  ,   gl  (ML) 

WHEN  NL  =  1 

<^  NL 

>  1 

YES 


CALCULATE 

PI  (NL)  ,   Ql (NL) 

WHEN  NL  >  1 


YES 


NO 


73 


0 


#25 


#26 


#27 


CALCULATE 
ARRAYS  AF(J) 
Ql(J) 


CALCULATE 
ARRAYS  BF(J) 
■AND  Q2 (J) 


CB(1)  ,   CB(2)  , 
CI  (NL)  ,  C2  (NR) 
ARE  CALCULATED 
IN  SUBROUTINE 
CENTER,  SEE 
SECTION  7.3 


YES 


#28 


CALCULATE 
C2(J)    J  =  NR 
+  1,   TO  J  = 
N2  -  1  AT 
TTIMPl 


74 


#29 


#30 


YES 


1 

NO 

r 

CALCULATE  Cl (J> 

J  =  1  tc  J  = 

ML  -  1 

AT  TIM£ 

ITIlMPl 

r 

YES 

f 

PRINT  Cl (J) , 

C2 (J) ,  CB (1) , 

CB (2) ,    KL,  NR, 

MM,  DTS 

,  TTMMP]. 

r 

NO 


75 


#31 


#32 


#33 


PRINT  AND  PLOT 
CONCENTRATIONS 
AT  PRESENT  VAL- 
UE OF  TTMMPl 


8 

<^     MM.  > 

KME 

YES 

r 

DIAGNOSTIC 

INCREASE 

WME 

UPVALUE  MM  AND 
PERFORM  NEXT 
ITERATION  OF 
i'lAIN  LOOP  78 


76 


0 


#35 


r 

PRINTOUT 
ARRAY 
PARAM  (J) 

#36 


#37 


GO  TO  73 

X 
©— 

f  END  j 


77 


8.2.     Flow  Diagram  for  Subroutine  CENTER 
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8.3.     Flow  Diagram  for  Subroutine  TRIDNL. 
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9.     HOW  TO  USE  DISTRIB 


9.1.  A  Brief  Description  of  DISTRIB:     VERSION  1 

The  intent  of  this  section  is  to  give  a  quick  general  introduction  to 
DISTRIB:     VERSION  1. 

The  meaning  of  the  undefined  terms  in  this  chapter  is  explained  in  the 
glossary,   chapter  12.     A  more  detailed  description  of  the  physical  as- 
pects of  the  redistribution  problem  is  given  in  chapter  2.     Chapter  3 
contains  a  more  detailed  description  of  the  overall  computational  pro- 
cedure. 

DISTRIB:     VERSION  1  calculates  a  redistributed  initial  impurity  profile 
in  a  crystal  wafer  that  has  undergone  an  oxidation  at  a  constant  tem- 
perature for  an  input  time  t  =  TFINAL.     The  impurity  concentration  in 
the  oxide  Ci   (2,t)   is  calculated  at  time  t  =  TFINAL.     Here,  z  denotes 

(fig.  3)   the  distance  measured  into  the  oxide  from  the  interface  be- 
tween the  ambient  oxygen  and  the  oxide  at  time  t  =  TFINAL.     (The  z- 
coordinate  system  changes  with  time.)     Also  calculated  at  time  =  TFINAL 
are  the  impurity  concentrations  in  the  silicon  C2(y/  t)  where  y  repre- 
sents  (fig.   3)   distance  into  the  oxide-silicon  composite  measured  from 
the  position  of  the  oxygen-silicon  interface  at  time  t  =  0.     At  this 
initial  time,  the  silicon  is  completely  unoxidized.     The  principal  out- 
put quantities  are  the  concentrations  CI  (L)   =  C^ ( (L  -  1)DZ,  TFINAL)  for 
L  =  1,  NL  where  NL  represents  the  index  of  the  first  grid  point  to 

the  left  of  the  position  of  the  moving  oxide-silicon  interfacial  boun- 
dary,  Zo(t),   at  t  =  TFINAL   (fig.   3).     Also  computed  are  C2(K)   =  C2 ( (K  -  1) 
DY,  TFINAL)   for  K  =  NR  +  1 ,    . . . ,  N2  where  NR  is  the  index  of  the  first 
grid  point  to  the  right  of  the  moving  boundary,  yo(t) ,  at  time  t  -  TFINAL 

(fig.   3)   and  N2  represents  the  number  of  grid  points  in  the  original  sili- 
con wafer.     The  quantities  DZ  and  DY  above  are  uniform  grid  widths  in  the 
z-  and  y-coordinate  systems.     Besides  the  concentrations  CI (  )   and  C2 (  ), 
DISTRIB  computes  CB(1)   =  Cj (Zq (TFINAL) )   and  CB(2)   =  C2 (Yq (TFINAL) ) ,  the 
impurity  concentrations  in  the  oxide  and  silicon  at  the  moving  oxide- 
silicon  interface  at  time  TFINAL.     All  of  these  concentrations  are  also 
plotted . 

9.2.  Definitions  and  Other  Information  Related  to  Input  and  Output 

In  section  9.21  the  real  input  quantities  are  briefly  defined.     A  table 
follows  the  definitions  giving  further  information  about  the  previously 
defined  input,   such  as:     the  computational  block  number  in  the  listing 
that  contains  the  read  and  format  cards,  format  information,  data 
sources,  and  numerical  ranges  for  which  the  parameters  have  been  suc- 
cessfully tested,  etc.     This  same  procedure  is  used  to  describe  the  in- 
teger input  in  section  9.22.     The  same  procedure  is  also  used  to  describe 
the  real  and  integer  output  in  sections  9.23  and  9.24.     Consult  the 
glossary  for  information  concerning  undefined  quantities  and  places  in 
the  documentation  where  more  thorough  definitions  are  given.     All  physi- 
cal quantities  must  be  entered  in  units  of  hours  and  micrometers. 
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9.2.1.     Real  Input 


The  real  input  is  defined  in  table  16  in  the  same  order  they  are  read 
(see  computational  block  #2  in  the  listing,  chapter  11,  for  the  read 
and  format  card) . 


Table  16:     Brief  Definitions  of  the  Real  Input 


Card  1 

CBULK      The  program  handles  automatically  piece-wise  constant  initial 
concentrations  in  the  silicon.     The  grid  points  in  the  silicon 
beyond  the  input  grid  point  NR  are  assigned  the  value  CBULK 
representing  the  concentration  in  the  bulk  of  the  silicon. 

CMAX        The  grid  points  1  through  NH  are  assigned  initial  concentrations 
of  magnitude  CMAX. 

ALPHA      A  given  voliime  of  silicon  swells  to  a  volume   (ALPHA)  times 

(the  original  amount  of  silicon  at  the  oxide-oxygen  interface) . 

TFINAL    The  total  time  that  the  wafer  is  to  undergo  redistribution. 

SEG  The  redistribution  coefficient. 


Card  2 

DFl  The  diffusion  coefficient  of  boron  in  the  oxide. 

DF2  The  diffusion  coefficient  of  boron  in  the  silicon. 

WFRC        The  fraction  of  a  mesh  width  distance,  DZ ,  that  the  moving 

boundary  moves  in  the  first  iteration  of  loop  78,  i.e.,  in  the 
first  time  step. 

BDRY        A  positive  number  representing  the  width  of  the  wafer. 

CONPRO    A  proportionality  constant  representing  the  evaporation  rate  of 
boron  from  the  oxide  into  the  oxygen  ambient. 


Card  3 

CONOUT    A  constant  representing  the  concentration  of  boron  in  the  oxygen 
at  the  oxygen-oxide  interface.      (In  our  testing  CONOUT  was  al- 
ways set  equal  to  zero.) 


CNSTA 1 
CNSTB [ 
CNTAU ) 


Constants  appearing  in  the  defining  equation  of  the  moving  boun- 
dary,  section  4.2. 


YIS 


Same  definition  as  Y2S,   etc.,  below. 


Card  4 

Y2S  I  Quantities  that  are  used  to  scale  the  concentrations  that  are 
X2S >  outputted  in  the  plotting  routines,  section  4.4,  and  computa- 
XIS)        tional  block  #34. 
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Runs  of  DISTRIB  have  been  successfully  completed  using  data  in  the  ranges 
shown  in  table  17.     Sources  for  the  data  and  formal  information  are  also 
given. 


Table  17:     Real  Input   (4  F0RI4AT   (5  F  10.5))    (see  block  #2  in  listing, 

chapter  11) .     Physical  quantities  are  assumed  to  be  in  units 


of 

micrometers . 

VARIABLE 

RANGE    (Note  4) 

(1) 

CBULK 

1.      (Note  4) 

(2) 

CMAX 

1.      (Not  tested  with  other  values) 

(3) 

ALPHA 

0.44 

(4) 

TFINAL 

0.1  to  4. 

(5) 

SEG 

(6) 

DFl 

Qop>  [41     rfii     r7i  [Ri 

(7) 

DF2 

(8) 

^'JFRC 

\J  *  \J  J-     L-SJ     W.JO  ^i>J^^Uti 

(9) 

BDRY 

o  .     uo    o  • 

(10) 

CONPRO 

n  1  t-r-i  1  non 

(11) 

CONOUT 

All    vnnQ    Vi^^7(^   }~if^^n    wi  t'h    'Zf^T'O  Vif^f^^ncio 

 L     J-  Lillo      LLCLVKIZ     XJ^dl     WXl^ii     ^K^d-\J  JOC^^QLiow 

of  our  experimental  conditions. 

(12) 

CNSTA 

See  tables  in  [1] 

(13) 

CNSTB 

^  :  

(14) 

CNTAU 

II 

(15) 

■   1 

YIS 

5. 

(16) 

Y2S 

2. 

(17) 

 _  1 

X2S 



4. 

(18) 

■ 

XIS 

2. 

Notes :      (1)     The  ranges  in  the  subsections  are  values  from  successful  runs 
and  are  not  to  be  interpreted  as  precise  limitations. 

(2)  This  relationship  is  cited  in  [4] . 

(3)  WFRC  must  be  chosen  small  enough  so  that  DM  =  WRRC*DZ  < 
ZCRIT.     For  an  explanation  of  this  restriction,  see  the 
comment  at  the  end  of  block  #7  of  the  main  program  in  chap- 
ter 7. 

(4)  The  problem  with  CONOUT  =0.   is  homogeneous  in  the  concen- 
trations and  therefore  all  runs  have  been  made  with  CBULKl. 
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9.2.2.     Integer  Input 

Table  18  contains  brief  definitions  of  DISTRIB's  integer  input. 

Table  18:     Brief  Definition  of  Integer  Input 

Card  5 

NZO  A  parameter  used  to  determine  the  magnitude  of  DZ,  eq  (4.3). 

MME  The  concentrations  are  calculated  at  times  ti,  tjyuyu^.  MME 

must  be  large  enough  for  tjvjME  >  TFINAL  or  the  program  will  stop 
before  calculating  the  concentrations  at  TFINAL.      (See  table  and 
sec.   9.3  for  guidance  in  choice  of  MME.) 

ND  This  integer  determines  how  many  grid  points  in  the  silicon  ad- 

jacent to  the  moving  boundary  will  be  plotted. 

NYO  The  number  of  mesh  widths  of  length  DY  that  equal  BDRY. 

NH  The  grid  point  on  the  y-axis  at  which  the  initial  boron  concen- 

tration changes  from  CMAX  to  CBULK,  eq  (2.10)  and  section  7.2, 
block  #10. 

KPRINT    The  concentrations  and  a  few  other  integer  parameters  such  as 

MM,  NL,  NR  will  be  printed  out,  whenever  MM  =  0  modulo  (KPRINT). 


The  range  of  values  of  integer  data  successfully  calculated  in  DISTRIB 
is  given  in  table  19. 


Table  19:     Integer  Input  Information 
Integer  Input   (1  FORMAT   (616))    (see  block  #2  of  listing) 


VARIABLE 

RMvIGE 

(1)  NZO 

2.   to  8. 

(2)  MME 

up  to  4000 

(3)  ND 

20.   to  30. 

(4)  NYO 

100  to  800 

(5)  NH 

10.   to  20. 

Certain  constants  that  can  be  used  to  control  various  aspects  of  how 
DISTRIB  runs  but  do  not  have  to  be  changed  from  run  to  run   (and  conse- 
quently, are  not  entered  formally)   are  listed  in  the  next  table. 
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Table  20:     Constants  Used  in  Program  But  Not  Made  Part  of  Input. 


(1)  0.2  and  100  used  in  block  #16  in  the  listing  of  the  main  program 
(see  sec.   7.1.,  block  #16). 

(2)  NMQD  and  NMODl   (sec.  4.4.)    and  block  #34  of  listing  of  the  main 
program.     These  quantities  need  to  be  changed  frequently  to  get 
plotted  points  nicely  spaced.     We  recommend  starting  with  NMOD  -  8, 
NMODl  =  4. 

(3)  BLTH   (sec.  4.3).     The  value  0.02  micrometers  has  worked  all  the 
time . 

(4)  In  printing  out  CI,  C2  in  blocks  #30  and  #34,  certain  indices  are 
skipped.     Of  course,  this  is  something  one  changes  from  run  to  run 
routinely,  in  order  to  find  the  values  at  the  grid  points  that  in- 
terests one. 


9.2.3.     The  Real  Output  of  DISTRIB 


Table  21:     Brief  Definitions  of  the  Real  Output.     See  Tables  24  and  25. 


DISTRIB  calculates  concentrations  in  the  oxide  and  silicon  for 
a  sequence  of  iterations  of  loop  78  with  indices  MI^l  =  1, 
I'lME .     The  time  at  the  end  of  the  MMth  iteration  of  this  loop 
is  TTi-UlPl   (in  units  of  BTM)  . 

The  concentrations  of  boron  in  the  oxide  and  silicon  at  the 
moving  boundary . 

These  are  the  concentrations  in  the  oxide  and  silicon  at  the 
Jth  grid  point   (for  the  specific  values  of  J  see  below)  at 
time  TTMMPl. 

Cl(l)         The  concentration  at  the  oxide-oxygen  interface. 
See  table  25. 

PARAM  is  an  array  that  contains  material  constants  (diffusion 
constants,  etc.)   in  the  original  units  of  micrometers  and 
hours.     SEG  is  the  redistribution  (=  segregation)  coefficient. 
The  swelling  constant  ALPHA  represents  the  reciprocal  of  the 
factor  that  a  given  volume  of  silicon  swells  when  converted  to 
oxide  at  the  moving  boundary.     TFINAL  is  the  length  of  time  of 
oxidation.     ZOCR  is  the  position  of  the  moving  boundary  in  the 
z-coordinate  frame  in  micrometers.     SCALE2  is  a  scale  factor 
used  in  the  plots. 


Further  specific  information  about  the  previously  defined  output  is 
given  in  table  22. 


PARAM 

SEG 

ALPHA 

TFINAL 

ZOCR 

SCALE 2 
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Table  22:     Real  Output   (5  FORMAT   (7F14.8))    (see  blocks  #30  and  #34). 


Whenever  MM  =  1  modulo   (KPRINT) ,  we  print  out 

(1)  C1(J)    for  J  =  1,  NL,  15. 

(2)  C2(J)    for  J  =  NR,  N2 ,  20. 

(3)  TTMMPl   (in  units  of  BTM) . 

(4)  Cl(l).   CB(1),  CB(2),  and  DTS   (in  BTM  units). 
In  addition,  when  TIMEND  <  TTMMPl,  we  print  out 

(5)  Array  PARAM,  SEG,  ALPHA,  TFINAL,  ZOCR,  and  SCALE2    (all  in  mi- 
crometers and  hours) . 

9.2.4.     The  Integer  Output 

Table  23:     Definitions  of  Integer  Output.     See  tables  24  and  25. 

NL  is  the  index  of  the  grid  point  in  the  oxide  that  is  nearest 
NL,  NR      the  moving  boundary  grid  point.     NR  is  the  grid  point  in  the 

silicon  that  is  nearest  the  grid  point  at  the  moving  boundary. 


MM 


The  index  of  loop  78  that  calculates  the  solution  at  time  t 
TTMMPl  from  the  known  solution  at  time  t  =  TTMM. 

See  table  25. 


NZO 
MME 
ND 
NYO 
NH 
KPRINTJ 


These  quantities  were  defined  in  table  11 


9.3.     A  Simple  Example 

In  this  section  we  illustrate  how  DISTRIB:     VERSION  1  is  used  to  com- 
pute the  redistribution  of  an  initially  uniform  distribution  of  boron 
under  wet  oxidation  conditions .     The  material  constants  used  in  this 
sample  calculation  were  chosen  to  obtain  the  best  fit  to  the  experimen- 
tal data  shown  in  figure  13.     For  further  comments  about  the  proper  val- 
ues for  these  material  constants,   see  the  end  of  this  section. 

The  first  card  in  computational  block  #2  reads  the  real  input  under  for- 
mat four.      (This  format  card  is  also  in  computational  block  #2.)  Since 
the  initial  boron  distribution  is  uniform,  the  value  of  CMAX  which  is 
the  concentration  at  the  first  NH  grid  points  in  the  silicon  (y-region) 
is  equal  to  CBULK  which  represents  the  initial  concentration  in  the 
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bulk  of  the  silicon  and  is  the  value  assigned  in  the  program  to  the 
grid  points  NH  +  1,  NH  +  2,   . . . ,  N2 ,     Because  the  partial  differential 
equations  are  homogeneous  in  the  concentrations  when  CONOUT  (sec.  2.1) 
is  zero,  as  in  this  example,  we  normalize  the  input  quantities  CMAX, 
CBULK  by  the  value  of  CBULK  (atoms  per  cubic  micrometer);  therefore, 
these  two  input  quantities  have  the  value  1  and  the  output  concentra- 
tions of  these  quantities  should  be  multiplied  through  by  CBULK  to  ob- 
tain the  true  concentrations.     Once  again,  the  glossary  is  for  answering 
these  questions.     The  next  parameter  read  is  the  swelling  ratio  ALPHA 
(sec.  2.1)   and  it  is  assigned  the  value  0.44  [4]. 

We  desire  to  find  the  concentration  of  the  redistributed  profile  after 
TFINAL  =  0.3  h.     When  TTMMPl   (T^mpi)    is  greater  than  TFINAL,   the  program's 
control  is  shifted  to  the  output  printing  blocks  #34  and  #35.     The  segre- 
gation coefficient  SEG   (dimensionless)   is  chosen  to  be  0.2,  in  order  to 
obtain  the  fit  shown  in  figure  13.     The  diffusion  coefficients  DFl,  DF2 
for  the  boron  in  the  oxide  and  silicon,  respectively,  were  specifically 
assigned  the  values  0.001  and  0.13  3  in  order  to  obtain  the  fit  shown  in 
figure  13,  and  these  values  are  within  the  range  of  the  reported  values 
for  these  quantities   (see  the  end  of  this  section) .     (Note  that  the  ra- 
tios of  the  diffusion  constants  used  in  the  program,  see   [7] ,  were  re- 
laxed in  order  to  obtain  the  fit  in  this  example.)     The  meaning  of  the 
real  grid  input  parameter  WFRC  =  0.5  is  explained  in  section  7.1  in  the 
material  concerned  with  block  #7. 

The  meaning  of  the  next  parameters  read,  CONPRO   (micrometers  per  hour) 
and  CONOUT   (atoms  per  cubic  micrometer) ,  are  related  to  the  flux  of  boron 
between  the  ambient  and  oxide   (sec.  2.1,  eq  (2.11)).     As  mentioned  al- 
ready, CONOUT  =  0  in  this  example  and  CONPRO  is  chosen  to  have  the  value 
1;  other  calculations  show  that  results  are  not  sensitive  to  the  choice 
of  CONPRO  including  the  value  zero.     The  next  three  parameters  read, 
CNSTA   (micrometers) ,  CNSTB   (square  micrometers) ,  CNTAU   (hours) ,  relate 
to  the  moving  boundary  motion  (sec.  4.2,  eq  (4.5)) .     The  wet  oxidation 
is  assumed  to  take  place  at  1100°C  and  therefore  from  table  1  in  [1] 
CNTAU(t)   =  0.0,  CNSTA(A)   =  0.11,  and  CNSTB(B)   =  0.51.     The  real  input 
read  next,  YIS ,  Y2S,  X2S,  and  XlS   (sec.  4.4.),  are  parameters  that  as- 
sign maximum  abscissa  and  ordinates  that  are  plotted  in  regions  one  and 
two,   i.e.,  the  oxide  and  silicon  regions,  respectively,  in  the  graphical 
output  in  computational  block  #34.     The  values  YIS  =  5,  Y2S  =  2,  X2S  =  4, 
and  XlS  =  2  have  always  yielded  good  graphical  output,  and  we  recommend 
their  use. 

The  integer  input  for  DISTRIB  is  read  following  the  real  input  in  com- 
putational block  #2  under  format  number  1   (also  in  computational  block 
#2).     The  integer  NZO  =  3,  whose  significance  is  explained  in  section 
4.1  below,  eq   (4.2),   is  read  first.     The  quantity  MME  =  1500  represents 
the  number  of  iterations  we  a  priori  assign  to  the  loop  78  (computational 
block  #14).     In  this  example,  only  595  iterations  of  loop  78  occur  before 
"^I^MPl  ^  TFINAL  =  0.3  h  and  the  program  stops  after  printing  the  output 
in  computational  blocks  #34  and  #35.     Integer  ND  represents  the  number 
of  grid  points  in  silicon  to  the  right  of  the  moving  boundary  point  in 
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the  silicon  that  are  plotted  in  computational  block  #34   (see  DO  LOOP  234) 
The  integer  NYO  -  350  is  read  next  and  its  function  in  defining  DY(Ay) 
is  explained  by  eq   (4.1).     The  integer  NH  =  12   (sec.   2.1,  eq  (2.13)) 
signifies  that  the  initial  concentration  value  CMAX  is  assigned  to  the 
first  twelve  grid  points  in  the  silicon  region  at  time  zero.     Since  in 
this  example  CMAX  =  CBULK,  it  makes  no  difference  what  the  value  of  NH 
is  except  NH  should  be  less  than  NYO  in  order  for  the  problem  to  make 
sense.     The  print-out  parameter  KPRINT  =  1000  signals   (in  block  #30)  the 
program  to  print  out  CI (  ) ,  C2 (  ),  the  concentrations  in  the  oxide  and 
the  silicon,   respectively,   at  the  end  of  the  iterations    (1  +  L*KPRINT) 
of  loop  78  where  L  =  0,  1,  etc.     Since  the  largest  value  of  MM  is, 

for  this  example,   595,  the  only  values  of  arrays  CI (  ),  C2 (  )  printed 
out  are  those  at  the  end  of  the  first  iteration  of  loop  78. 

These  are  certain  constant  parameters  that  are  not  read  into  the  program. 
These  are  described  in  section  9.2   (table  20). 

For  various  literature  references  for  material  constants ,  we  recommend 
[1,3,4,5,5,7,9,12,13]   in  the  supplemental  bibliography  in  chapter  13. 

The  output  of  DISTRIB  is  listed  and  described  in  section  9.2.     The  print- 
out in  block  #30  is  effected  after  every  MM  =   (1  +  L'KPRINT)th  iteration 
(for  L=0,  1,   2,    ...)  of  loop  78,  which  is  the  loop  computing  the  con- 
centrations at  the  successive  times  t^M.     Since  KPRINT  -  1000  for  this 
example  and  only  595  iterations  of  loop  78  are  required  to  compute  the 
concentrations  at  time  =  TFINAL  =  0.3  h  in  this  example,   the  only 

concentrations  printed,  under  block  #30,  are  those  after  the  first  itera- 
tion and  are  shown  in  table  24.     From  this  table  and  computational  block 
#30,  we  see  printed  out  the  values  of 

CI (J) ,  J  =  1,  NL,  15 

C2  (J)  ,   J  =  NR,  N2,  20. 

Since,  at  time  t^^  =  0.0389   (BTM  units)  ,   there  is  only  one  fixed  grid 
point  in  the  oxide,  the  only  value  in  the  oxide  that  is  printed  is  the 
grid  point  at  the  oxide-oxygen  interface.     However,  between  the  values 
of  Cl (  )   and  C2 (  )   that  we  have  printed  out  in  table  24,  there  appear 
the  values  of  boron  concentration  at  the  moving  boundary  grid  points  in 
the  oxide  and  silicon  CB(1),  CB(2).     Also  printed  in  the  same  row  with 
these  quantities  are  DTS   (in  BTM  units)   and  the  mesh  width  lengths  in 
the  oxide  and  silicon,   respectively,  DZ  and  DY  in  units  of  BLTH  =  0.02 
ym.     In  the  final  print-out   (table  26),  BTM  =  0.0030  h  is  printed.  Also, 
table  24  shows  the  last  index  MM  of  loop  78  and  index  NL  of  the  first 
fixed  grid  point  to  the  left  of  the  moving  boundary  in  the  oxide.  Final- 
ly,  in  the  same  line  appears  the  index  NR  of  the  first  grid  point  to  the 
right  of  the  moving  boundary  in  the  silicon. 

The  printout  from  block  #34  is  shown  in  table  25.     The  main  output  are 
the  concentrations 

Cl  (J) ,  J  =  1,  NL,  15 

C2  (J)  ,  J  =  NR,   N2,   20.  :. 
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These  are  the  concentrations  at  the  first  time  that  TTMMPl   (=  99.756). 
is  greater  than  TFINAL   (in  BTM  imits) .     It  is  seen  from  table  25  that 
this  occurs  at  the  end  of  the  MMth  =  595  iteration.     The  meaning  of  the 
other  output  in  table  25  is  similar  to  the  corresponding  output  in  table 

24  that  has  already  been  described   (see  also  sec.  9.23). 

The  concentrations  at  time  TTMMPl  =  99.756   (in  BTM  vmits)   or  TTMMPl  ~ 
TFINAL  =  0.3  h  in  the  oxide  are  plotted  in  figure  11.     The  distances 
from  the  oxide-oxygen  interface   (scaled  in  ZOPl  units  which  is  the 
length  of  the  oxide  region  at  time  TTMMPl  ~  TFINAL)   are  plotted  on  the 
X-axis.     The  concentrations  at  the  positions  of  the  first  and  (1  + 
L*NMOD)th  grid  points   (for  L  =  1,   2,    ...)   are  plotted  on  the  y-axis,  in 
this  example,  NMOD  =  8.     The  value  of  ZOPl   (in  micrometers)   is  computed 
in  the  last  table   (26) .     At  the  same  time  that  the  concentrations  in  the 
oxide  are  plotted  in  figiare  11,  the  concentrations  in  the  silicon  are 
plotted  in  figure  12.     The  x-axis  in  this  figure  represents  distance  in 
the  unit  UN2  =  2SQRT  (DF2 'TTI-^MPl)    (where  TTMMPl  is  very  nearly  equal  to 
TFINAL)   from  the  oxide-silicon  moving  boundary.     The  concentrations  at 
the  moving  boundary  (in  the  silicon)   and  at  the  location  of  the  grid 
point  with  index  NR  are  plotted.     After  this,   the  concentrations  at  the 
grid  points  with  the  successive  locations    (NR  -  1)   DZ  +   (J  -  2) 'NMODl'DZ, 
for  J  =  3,    .  .  .  ,  ND   (with  NMODl  =^  4  in  this  particular  example)  are 
plotted.     The  lanit  length  of  the  abscissa,  UN2   (micrometers),  is  com- 
puted in  table  26  from  computational  block  #35.     In  addition,   in  table 
26  are  shown  all  the  input  quantities  in  units  of  micrometers  and  hours 
as  printed  out  in  computational  block.  #35.  Besides  the  input  data  that 
are  printed  out,   the  quantities  ZO (TFINAL)    (this  is  the  same  as  ZOPl, 
when  TFINAL  =  TTMMPl)   and  2 'SQRT (DF2 'TFINAL)    in  micrometers  which,  as 
we  have  already  explained,  are  used  as  quantities  to  scale  the  x-axis 
in  the  oxide  and  silicon  concentration  plots  in  figures  11  and  12.  Also 
printed  is  the  value  of  BTM  in  hours;  hence,   the  times  in  tables  24  and 

25  which  are  in  units  of  BTM,  can  be  converted  to  hours,   if  so  desired. 
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Table  24.     Example  of  output  of  concentration  in  oxide  and  silicon  at  end  of 
iteration  MM  =    (1  +  L*KPRINT)    in  loop  78  for  L  =  0. 


C  I  <  .  )    AT    T IME 

2. 52280046 
NL 
1  1 
CI  (  I  ) 
2  .522eC046 
)    AT  TIME 
.  9a4e'*473 
1  .OOCCCOOG 

1  .oocooooo 


MM 


C2(  , 


.038^98095 


CB  (  I  ) 

3. 68795223 
.033996095 

1 . 00  cooooo 
I. cooooooo 

1.00000000 


CB(  2) 

.77759045 

1  .00000000 
1  .00000000 

1. ocoooooo 


DTS 

. 03899809 

1 .00000000 
1 . 00000000 
I .00000000 


DZ 

.541 12554 


1  .00000000 
1  .00000000 


DV 

.71428571 


1  .00000000 
1 .00000000 


I . 00000000 
1 . 00000000 


Table  25.     Example  concentrations  in  oxide  and  silicon  when  TTMMPl  =  TFINAL. 


C  l(  .  )    AT    T  IME 

.03302989 
UH  NL 
595  32 
C  1  (  1  ) 
.03302989 
C2<  .  )    AT  TIME 
.  189  53686 
1 .00000000 
1 .00000000 


99.756364064 
.8754271 7 


NR 
12 


CB(  1  ) 

.813071 12 
99.756364064 

. 844  37386 
1 .00000000 
1. OCOOOOOO 


, 83500469 


CB(  2) 

.16261422 

. 98733067 
1 .00000000 
I  .00000000 


. 82 160333 


DTS 

. 28899809 


. 99958049 
1  .00000000 


.81678325 


OZ 

,541 12554 


.99999438 
1  .00000000 


■81453516 


DY 

.71428571 


.99999997 
1. OOOOOOOO 


.81331455 


1  .00000000 
1 . OOOOOOOO 


Table  26. 


ITA    FOH  STEP 
OF  1 
. 00 1 0000  0 

SEG 
.20000000 
CB<  1  ) 

.81 3071 1 2 
CMAX 
1  .OOOOOOOO 

Nzc  fue 

3  1500 


Listing  of  input  physical  data  (micrometers, 
computational  data. 


liours)   and  important 


NO 
30 


1 

DF2 
.  1  33000  00 
ALPHA 
.44000000 
CB(  2  ) 
.  1626  1422 
CBULK 
1.00000000 
NYO 
350 


CONPRQ 
1  .OOOOOOOO 
WFRC 

.05000000 
FI RST+DT 

. 0008691 7 
VIS 

5.00000000 
KPR I  NT 
1000 


CONOUT 
.OOOOOOOO 

TFINAL 
. 30000000 
REAL  TIME 
. 30001914 

y2S 

2.00000000 


CN5TA 
.  I  1 000000 
ZO ( TF I NAL ) 
. 34001 236 
BDRV 
5  .OOOOOOOO 

X2S 
4  .OOOOOOOO 


CNSTB 
. 5 1 000000 
2*SQRT{ DF2*TF INAL ) 
.39949969 


CNTAU 
, OOOOOOOO 


XI  S 
2. OOOOOOOO 


BTM 

.00300752 
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Figure  13 .     Comparison  of  three  calculated  impurity  concentration 
distributions  in  silicon  with  experimentally  determined  concentra- 
tion data. 
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10.     VERIFICATION  OF  PROGRAM 


10.1.  Validation  of  the  Program 

A  solution  calculated  by  DISTRIB  for  a  problem  with  the  specified  con- 
centration C2^(0,t)   =  Cq  at  the  oxygen-oxide  interface  has  been  compared 
with  an  analytic  similarity  solution  given  in   [4] .     The  calculated  boun- 
dary concentrations  C2^(zQ(t),  t)   agree  with  the  similarity  solution  [5] 
to  within  a  few  percent  for  a  variety  of  problems  with  differing  segrega- 
tion constants.     The  calculated  solutions  also  have  qualitative  proper- 
ties possessed  by  the  similarity  solution  such  as  moving  boundary  con- 
centrations that  are  constant  with  time.     There  is  also  an  analytic  solu- 
tion to  special  problems  with  zero  flux  at  the  oxygen-oxide  interface. 
In  this  problem,  the  concentration  in  the  oxide  or  z-region  is  constant 
and  DISTRIB 's  calculated  solution  is  constant  also  for  these  conditions. 

10.2.  Limitations  of  the  Program 

The  previous  remarks  must  all  be  qualified  for  very  short  times,  up  to 
5At   [5] .     The  difficulties  in  obtaining  the  solutions  for  the  first  few 
iterations  are  similar  to  those  described  in   [2]   for  analogous  problems. 
The  problem  with  the  inaccuracy  of  the  first  few  iterations  can  be  some- 
what mollified,  however,  by  choosing  smaller  and  smaller  values  of  WFRC. 
Our  recommendation  for  obtaining  a  solution  at  small  times  is  to  calcu- 
late solutions  using  a  series  of  decreasing  values  of  WFRC  and  then  use 
those  solutions  to  extrapolate  the  solutions  to  time  zero. 
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11.      PROGRAM  LISTINGS 


1.     DISTRIB:     VERSION  1  MAIN  PROGRAM 


1 

5 

41 
1+3 

8 
9 

26 
27 
36 

32 
39 
4b 


*+0 
38 

28 

29 
30 
31 
33 

35 

20 

21 
22 
23 

24 
25 
19 
lb 


START  COVPUTATIOM  BLOCK  f^UMBEP  1 
IMPLICIT  POUBLE  PRECISION  (A-H»0-Z) 
REAL     XPL(200)  »YPL(200) 

DIMENSION  TT(2000)»ZO(20no)»AF(200n)rRF(2ono)»PARAM(21) 
DIMENSION     BT(20) »    DFA(2C) »nFB(?0) »CONPR(20) fCNST] {2Q)  »CMSt?(2o) » 
1   CNTA(2C)  .TFIfNlN(20)  »ZONEW(?n)  rTNEW(20) 
COMMOM/RL?0/DT»nY»DZ»THETA 
C0r^M0N/PL21/DFl  »DF2 
C0'^M0N/eL22/DZP  rPYP 

C0MMGN/BL23/  82M  »  RiiR  »  P2L  »  A2L  »  A2fw  ,  A2R  '  A I L  »  A  1  M »  ftlP»PlL»RlMrPlP 
C0MM0N/PL2'+/C2L»C2M»C2R»CIL»C1'^»C1R 

C0MM0N/BL25/P(2) »0(2)»PA(2)»QA(2)»C(2)fCP(2)»P(?)»S(2)»SM{P) 
C0MM0M/FL26/P1 (2000) »P2 (2000) fOl (2000) »Q2(?00n)»Dl (2000) 'D9(2nno) 
C0N'M0N/PL2B/N2M1  »NLM1  rM2    f  NRPl»NL»NR 
C0MM0rVBL29/Cl (2000) ►C2 (2000) 

COMMON  /BL30   /  CO  f  CK  »  TZERO  »  JEE »  ALPHA  f  C  AP  A  »  SeG  »  CWY 

COMMON  /BL31   /  PIST(100)»    SOL (100) 

C0MM0N/PL32/C0MPRO»C0N0UT 

END         COMPUTATION  BLOCK  ^lUMBER  1 

START     COVPUTATIOM~BLOCK   NUMBER        2  ~ 

READ ( 5r4) CRULK  »CMAX» ALPHA »TFlNALfSEG»DFl »DF2»wFPC»nDRY»C0NnR0» 
1C0N0UT»CNSTA»CNSTR»CNTAU»  Y1S»  Y2S»X2S»X.1S 
PEA0(5»  1  )NZ0»MME»MD»NY0»NH»KPRIr|T 
FORMAT  (616) 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

FORMAT  (16H  Cl(.)  AT  TIME  »1F14.9) 
FORMAT  (16H  C2(.)  AT  TIME  »1F14.9) 
FOPMAT(    IHl  ) 

FORMATCiOH  TLIN  CMTAU  TT^^'v  SLOPeZCRIT  CnSTX  ) 


(5F10.5) 
(7F14.8) 

(    4 OH  MM  ML  NR 

(40H  MME  MUST  BE  INCREASED 
( 16H  CI ( . )    AT  TIME 
)    AT  TIME 


FORMAT 
FORMAT 
FORMAT 
FORMAT 


 C/CB  VS  Z/ZO—  STEP  NO     flI5  ) 


(   45H  PROFILE   IN  OxIPE 
(IHl  ) 

(    40H  TIMENO  TNEW(J)    BTM     CNTAU  TLIN  BT(J)  ) 
(lOlH  CKl)  CP(1)  CR{2) 

DZ  OY 
COMPUTATION  BLOCK  NUMBER  2 

(    52H  Z0Pl»0F2»UN2»nYfnZfFNC»  (NL-1)  fTTMMpi,c;cALE2»XPSCUT 
(      19H  DATA  FOP  STEP  NO  »1I5) 
33H  NZO  MME  ND  NYO  NH  JSTOP 
43H  S(l)   ST  STM  El  E2  E3  ERROR 


) 


) 


) 
) 

CONPRO 


CONOUT 


DTS 
) 

END 

FORMAT 
FORMAT 
FORMAT( 
FORMAT ( 

FORMAT(     33H  LIST  OF  Z0( , ) 

FORMAT(      33H  LIST  OF  TT(.) 

FORMAT    (lOOH  DFl  DF2 

1  CMSTA  CNSTB  CNTAU  ) 

FORMAT{  67H  PROFILE  IN  SILICON  C/CP  VS  ( Y-YO ) / ( 2*SQRT ( DF2*TF INaL) 
1)     —  STEP  NO       »1I5  ) 

FORMAT    (21H  LIST     AF ( . )  ) 

FORMAT   (21H  LIST     PF(.)  ) 

FORMAT   (   30H  COXIDE  CBULKZOl   ALPHA  TFINAL  ) 

FORMAT  (103H  SEG  PRINT  DFl  DF2  WFRC  BDPY  CMA^  CONPRO  CONOUT 
1  CNSTA  CNSTB  CNTAU     YIS  Y2S  X2S  ) 

FORMAT (    33H  LIST  N2  NH     Nl   NR     NL  ) 

FORMAT(    33H  LIST  OZ  DY     OT  TT(1)   SM(1)    5M2)  ) 

FORMAT    (21H  LIST     PK.)    02(.)  > 

FORMAT   (   33H  LIST  tJ2L  B2P  P2M  PlL  RlR  ^IM  ) 
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17     FORMAT  (    33H  LIST  A2L  A2M  AlL  AIR  AIM  ) 

15     FORMAT(  20H  LIST  OF  PlOfQK.)  ) 

It^     FORMAT  (60H  VALUES  OF  RM  R7     RY  B  C  F  CMVl   CMV2  PIS  P?5  AT  m= 

1  »1I9  ) 

13     FORMAT  (61H  VALUES  OF  FNL  07P  DYP  TES  ZOPl   FN"     TT(MMPl)   SFS  AT 

1     MM=  fllB) 

12       FORMAT  (28H  DIVIDE  FAULT  AT  MOV.   RDRY.  ) 

11  FORMAT  (28H  DIVIDE  FAULT  AT  LOOP  107  ) 
10       FORMAT   (2PH  DIVIDE  FAULT  AT  LOOP  105  ) 

34       FORMAT   (lOlH  SEG  ALPHA  WFRC  TFlNAL 

1  ZO(TFINAL)        2*S0RT(nF2*TFINAL)  ) 

6       FORMAT   (   7F10.5  ) 

37       FORMAT(   80H         CBd)  CB(2)  FIRST+DT 

1  BDRY 

41       FORMAT   (lOlH  CMAX  CBULK  YIS 

1  X2S 
42     FORMAT    (  50H 


NH  KPPINT 


REAL  TIME 
) 

Y2S 
) 

) 


CMAX  CBULK 
XIS  BTM 
NZO       MME       ND  NYO 
C  DATA  FIXED   IN  VERSION  1 

DATA       DFACl) »DFP(1) »   CONPR ( 1 ) » CNSTl ( 1 ) » CNST2 ( D ^   CNTA ( 1 ) » TF IMN { 1) 
1     /.0001D0» .OOOIDO^    .OlDOf    . 027D0 » . IDO » . iDO » . 0 IDO/ 
TW=0, 
ZW=0. 
JSAVE=0 
JSTOP  =1 
ZOl=0. 
C0XIDE=1. 

C  JSTOP  TELLS  NUMBER  OF     DRIVE   IN  STEPS 

JSTM1=JST0P-1 

C  START     COMPUTATION  BLOCK  NUMBER  3 

5CALE2  =  2.*  SQRT(DF2*TFINAL) 
C  THESE  MATERIAL  PARAMETERS  ARE  RESCALED  THUS  WE  HAVE  TO  STO^E 

C  THESE  GIVEN  PARAMETERS   IN  UNITS  OF  MICRONS  ANP  HOURS 

PARAM(1)=DF1 

PARAM(2)=nF2 

PARAM{3)=  CONPRO 

PARAM(a)=CONOUT 

PARAM(5)=CNSTA 

PARAM(6)=CNSTB 

PARAM(7)=CNTAU 
C  END         COMPUTATION  BLOCK  NUMBER  3 

C  START     COMPUTATION  BLOCK  NUMBER  4 

BLTH=  .02 

C  DEC  MODIFIES  PROGRAM     TO  DIMENEnTIONALISE  ABOUT  DFC+DF2 

DFC=  1. 
DF2=DF2*DFC 
BTM=(RLTH**2)/DF2 
C0NPR0=C0NPR0*(BTM/BLTH) 
CNSTA=CNSTA/BLTH 
CNSTB=CNSTR/DF2 
CNTAU=CNTAU/BTM 
BDRY=RDRY*(1./BLTH) 
Z01=Z01* (1 ./BLTH) 
DF1=DF1/DF2 
DF2=1 ./DFC 

TIMEND=TFINAL*(1./BTM) 
C  END         COMPUTATION  BLOCK  NUMBER  4 

C  ZOl   IS  POSITION  ON  Z-AXIS  WHERE  THE  NON  DEGENERATE  {Z0{2ER0) 

C  NOT  EQUAL  TO  ZERO   )   MOVING  BOUNDARY  BEGINS 

C  ZCRIT   IS  The  position  on  the  Z-AXIIS  where  the  moving  BOUNDARY 

C  CHANGES  FROM  LINEAR  TO  PARABOLIC 

C  CALCULATE  DY»DZ»DT  AND  AUXILLARY  QUANTITIES 

Y01=ALPHA*Z01 
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C  START     COMPUTATION  BLOCK  MUI^BER  5 

DY=BDRY/NYO 

FNZO=FLOAT(NZO) 

DZ=DY*(1 ./( ALPHA*FNZO) ) 
C  END         COMPUTATION  BLOCK  NUMBER  5 

C  ================r==========r=r==r====r==========r=========r===r=== 

C  START     COMPUTATION  BLOCK  NUMBER  6 

CNA=,5*CN5TA 

cnb=cnstb 

c        take  out  following  card  if  cnstx  is  known 

CNSTX=SORT (CNA**2+CNB*CnTAU)  -CNA 
C  END         COMPUTATION  BLOCK  NUMBER  6 

CNTaU=CNSTX* (CNSTX+CNSTA)  /CNSTb 

C  START     COMPUTATION  BLOCK  NUMBER  7 

DM=DZ*WFRC 

C  INSERT   IF  MOVING  BDRY   IS  LINEAR   AND  THEN  PARABOLIC 

TLlNrCNTAU 

C       IF  MOVING  BOUNDARY  STARTS  AT  Y»0.   THESE  CARDS  TAKE  CARE  OF  STTUaTIO 
C  N  THAT  LINEAR  PART  OF  MOVING  BOUNDARY  IS  MOT  PPESFNT 

IF(   TLIN   .LE.    .00000001   ) DTr ( DM* ( DM+CNSTA )) /CNSTB 

IF(   TLIN   .LE.    .00000001  )nTS=DT 

IF{   TLIN   .LE.    .00000001  )ZCRIT=0. 

IF(   TLIN   .LE.    .00000001  )SL0PE=1. 

IF(TLIN   .LF.      .00000001   ) GO  TO  216 

ZCRITrCNSTX 

SLOPE=ZCRIT/CNTAU 

DT=DM/SLOPE 

DTS=DT 
216  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  7 

C  START     COMPUTATION  BLOCK  NUMBER  B 

PARAM( 10)=DT*BTM 

C  END         COMPUTATION  BLOCK  NUMBER  8 

C  START     COMPUTATION  BLOCK   NUMBER  Q 

C  CALCULATE   INITIAL     INITIAL  DISTRIBUTIONS 

DUM=Y01/DY 

NR=INT(DUM)+2 

N1=INT{Z01/DZ)+1 

NL=N1 

N2=NY0+1 

N2M1=N2-1 

C  END         COMPUTATION  BLOCK  NUMBER  9 

DO  98  J=1»N1 
CI ( J)=COXIDE 
98  CONTINUE 

CB{1)=C0XIDE 

C  ===r===============r==r===r=====r=r===r========r=======rr===r===r= 

C  START     COMPUTATION  BLOCK  NUMBER  10 

N1P1=M1+1 
NH=NH  +rJR 
NHP1=NH+1 
DO  96  J=NR  »NH 
C2( J)=CMAX 
96  CONTINUE 

DO  97  J=NHP1»N2 
C2( J)=CBULK 
97  CONTINUE 

CB(2)=CMAX 

C  END         COMPUTATION  BLOCK  NUMBER  10 
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START     COMPUTATION  BLOCK  MUMBER  11 
IF(Z01    .LF.     ZCRIT)   TTMM  =  701/5L0PE 


IF{Z01    .GT.     ZCRIT)   TTMM-     ( ZOl* ( ZOl+CMSTA ) ) /CNSTR 

END         COI^PUTATION  BLOCK  ^)UlV!BER 

11 

- ------ 

START     COMPUTATION  BLOCK  NUMBER 

^^37  —  — — —  —       r^  —  —          —  —  —  —  ^  —       —  ~  — 

FNL=FLOAT (ML) -1 . 

FNR=FL0AT(NR)-1. 

5M{1)=(Z01-FNL*DZ)/DZ 

SM(2)=(FNR*DY-Y01 ) /DY 

END         COMPUTATION  BLOCK  NUMBER 

12 

START     COMPUTATION  BLOCK  NUMBER 

13 

MAIN  LOOP 

P2(N2)=n. 

Q2(N2)=C2(rJ2) 

END         COMPUTATION  BLOCK  NUMBER 

13 

================================ 

START     COMPUTATION  BLOCK  NUMBER 

—  —  —  —  —  — —  —  —  —  ~  — —  ~  —  Z.  — —  ~  —  —  —  — —  —  —  ~  — —  —  —  —  — — 

ZOPS=0, 

ZOP1=0. 

TJ=0. 

N0RUN=1 

KPOW=0 

JSAVE=0 

DO  78  MM=1»MME 

END         COMPUTATION  BLOCK  NUMBER 

1*+ 

START     COMPUTATION  BLOCK  NUMBER 

15 

KW=MOD(MMrKPRINT) 

END         COMPUTATION  BLOCK  NUMBER 

15 

^  —  ZZ~^  —  ~  —  —  —  ~  —  —  ~~~~  —  ~  —  —  —  —  —  ~  — —  —  —  —  —  — 

START     COMPUTATION  BLOCK  NUMBER 

—  —  ———— ~~——~~— ~— —  —  —  ———————— ——————— — 

16 

SPEEDS  UP  MOVING  BOUNDARY 

ZOACEL=(ZOP1-ZOPS)/DZ 

IF(ZOACEL   .GE.      .2   )   GO  TO  88 

KP0W=KP0W+1 

IF(KPOW   .GT.    100  )DTS=DTS+,05 

IF(KPOW   .GT.   100  )DT=DTS 

IF(KPOW   .GT,   100   )CALL  TRTDNL 

IF(KPOW   .GT.    100   ) KPOW=KPOW-100 

CONTINUE 

ZOPS   IS  ZO(TTMM) 

ZOPS=ZOP1 

END         COMPUTATION  BLOCK  NUMBER 

16 

================================ 

START     COMPUTATION  BLOCK  NUMBER 

17 

IF(TJ   .GE.   1.)  NL=NL+1 

IF(TJ   .GE.  l.)Cl(NL)=CP(l) 

IF(TJ   .GE.   1. )FNL=FNL+1 . 

NLM1=NL-1 

NLM2=NL-2 

IF(TJ   .GE.  l.)DT=DTS 

IF(TJ   .GE.   1.)     CALL  TRIDNL 

IFCTJ   .GE.l.    )  SM(1)=0. 

END         COMPUTATION  BLOCK  NUMBER 

17 

START     COMPUTATION  BLOCK  NUMBER  18 

MMP1=MM+1 

TTMMP1=TTMM+DT 

IF(    TTMMPl      .LE.    TLIN   )      ZOPl   =    (    TTMMPl   -  TW)*  SLOPE  +ZW 
IF(TTMMP1    .GT.   TL IN ) Z0P1=SQRT ( ( CNA**2 ) +CNB* ( TTMMP l-TW ) )-CNfl+ZW 
Y0P1=ALPHA*Z0P1 

END         COMPUTATION  BLOCK   ^IUMBEP  18 
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C  START     COMPUTATION  BLOCK  NUMBER  19 

C  CHECK   IF  NEW  CELLS  ARE  ADDED 

SES= ( FNR+DY-YOPl) /DY 

IF(SES  .LT.   0.)  NR=NR+1 

IF(SES   .LT.   0.)  FNR=FNR+1, 

IF(SES   .LT.  0.)SM(2)=1. 

NRP1=NR+1 

NRP2=NR+2 

NRM1=NR-1 

C  END         COMPUTATION  BLOCK  NUMBER  19 

C  START     COMPUTATION  BLOCK  NUMBER  20 

MMP2=MMP1+1 
TTMMP2=TTMMP1+DT 

IF(   TTMMP2   ,LE.   TLIN   )   Z0P2=(TTMMP2  -TW   )*  SLOPE     +  ZW 

IF(   TTMMP2   ,6T,   TLIN   )   zOP2=SQRT { { CNA**2 ) +CNB* ( TTMMP2-TW ) ) -CNA+zW 

TJ= (Z0P2-FNL*DZ) /DZ 

IF(TJ   .LT.   1.)   GO  TO  219 

Z0P1=(FNL+1. )*DZ 

IF(   ZOPl   .LE.  ZCRIT)TTMMP1=(Z0P1-ZW)/5L0PE  +  TW 

IF(   ZOPl   .GT.   ZCRIT)TTMMP1=( (ZOPl-ZW)*(   ZOPl-ZW+CMSTA ) ) /CNSTB  +TW 

DT=TTMMP1-TTMM 

CALL  TRIDNL 

Y0P1=ALPHA*Z0P1 
219  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  20 

C  CONSERVATIOfJ  OF  NUMBER  CALCULATION 

S1M=C1  (rJL)*DZ*.5+(Cl(NL)+CB(l)  )*SM{1)*.5*DZ 

IF(NL   .EQ.   1)   S1M=(C1<1>+  CB ( 1 H *SM ( 1 ) * . 5*0Z 

IF(NL   .GE.  2)51rv'=SlM+.5*nz*Cl  (1) 

IF(NL   .LE.   2)   GO  TO  222 

DO  111  J=2»NLM1 

S1M=S1M+DZ*C1 ( J) 

111  COrJTINUE 

222     52M=C2(NR)    *DY* . 5+ ( C2 ( NR ) +CB ( 2 ) ) *  SM ( 2 ) * . 5*DY 
DO  112  J=NRP1»N2M1 
S2M=S?M+C2( J)*DY 

112  CONTINUE 

C  ======r==========r====r================rr=======r=====r===r==r==r= 

C  START     COMPUTATION  BLOCK  NUMBER  21 

5(1 )=(Z0P1-FNL*DZ)/DZ 
5{2)=(FNR*nY-Y0Pl )/DY 
C  END         COMPUTATION  BLOCK  NUMBER  21 

C  =r=======r=====r==z==r===========r======r===r==rr==r============== 

C  START     COMPUTATION  BLOCK  NUMBER  22 

IF(MM   .EQ.   1)   CALL  TRIDNL 
C  END         COMPUTATION  BLOCK  NUMBER  22 

C  START  COMPUTATION  BLOCK  NUMBER  23 

IF(NL.GT,    1)   GO  TO  223 

IF(NL  .EG.  1)DM1=S(1)*DZ 

IF(NL  .EQ.  l)DM2rDMl*DMl/(2.*DT) 

IF(NL  .EQ.   1    .AND. MM   .EG.     1)  DM2=0. 

IF(NL  .EG.  1)DM3=DM2 

IF(NL  .EG.   l)DM=l./(  DF1+0M2+DM1*C0NPR0) 

IFCNL  .eg.   1  )P1(1)=DM*DF1 

IF(NL  .EG.  1)Q1{1)=DM*(C1(1)*DM3+DM1*CONPRO*CONOUT) 
C           END         COMPUTATION  BLOCK  NUMBER  23 
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C  START     COMPUTATION  BLOCK  NUMBER  24 

223     IF(NL   .6T.   1    ) DM=1 . / {DZ*DZ/ (2 . *DT )   +DZ*  CONPRO+DFl  ) 
IF{NL  .GT.   1  )P1(1)=DF1*0M 

IF(NL   .GT.   1  )Q1(1)=(+C1(1)*(DZ**2)/(2.*Ot)+DZ*CONPRO*COMOUT1*Dm 
C       ^  END         COMPUTATION  BLOCK  NUMBER  24 

C  START     COMPUTATION  BLOCK  NUMBER  25 

IF(NL   .LE.  2)   GO  TO  213 
DO  203  J=2»NLM1 
AF( J)=A1M*C1( J) 

203  CONTINUE 

DO  205  J=2»NLM1 
JM1=J-1 

Q1U)  =  (AF(  J)-Q1(  JM1)*A2L)/D1  (  J) 
205  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  25 

IF(MM  .EQ.   KD)   WRITE(  6»15) 

IF(MM  .EQ.KD)  WR ITE ( 6 » 5) ( Pi ( J) » J=l f NLMl » 1 ) 
IF(MM  .EQ.KD) WRITE(6»5)  (Ql ( J) » J=l rNLMl » 1 ) 
IF(  MM  .EQ.  KD  )  WRITE(6,20) 

IF(MM  .EQ.  KD) WRITEC6»5)    (AF(J)  »J=2»NLMlfl) 
C  ======================r==============r======r========r==r===r===== 

C  START     COMPUTATION  BLOCK  NUMBER  26 

213  CONTINUE 

DO  204  d=NRPlfN2Ml 
BF( J)=B1M*C2( J) 

204  CONTINUE 
NE=N2-NRP1 

DO  207  J=1»NE 

K=N2-J 

KP1=K+1 

Q2(K)=(BF(K)-Q2(KP1) )/D2(K) 
207  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  26 

C  ==============================r====r========r===========r===rrr=r= 

C  START     COMPUTATION  BLOCK  NUMBER  27 

CALL  CENTER 
C  END         COMPUTATION  BLOCK  NUMBER  27 

C  START     COMPUTATION  BLOCK  NUMBER  28 

DO  209  J=NRP1»N2M1 
JM1=J-1 

C2(  J)=P2(  J)*C2(  JMl)-«-G2(  J) 
209  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  28 

C  START     COMPUTATION  BLOCK  NUMBER  29 

IF(NL  .LE.   1)   60  TO  215 
DO  210  K=1»NLM1 
J=NL-K 
JP1=J+1 

C1(J)=P1{J)*C1(JP1)+Q1(J) 
210  CONTINUE 
215  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  29 

C       ■    START     COMPUTATION  BLOCK  NUMBER  30 

IF(KW  .EQ.  0  .OR.     KW     .GT.   1)     GO  TO  799 
WRITE(6»32) 

IF(KW   .EQ.   1)  WRITE(6»8)  TTMMPl 

IF(KW     .EQ.   1     )WRITE(6»5)    (CKj)    »J=lrNL»l5  ) 
WRITE(6»41) 

IF(KW  .EQ.   1)  WRITE(6»1)  MM»NL»NR 
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WRITE(6»45) 

WRITE (6» 5)   Cl(l) »C6(1) »CR(2) »DT5»DZ»DY 
IF(KW   .EQ.   1)   WRITE{6»9)  TTMMPl 
IF(KW   .EQ.   1)WRITE(6»5)    (C2(J)»J=NR  »M2'?0) 
799  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  30 

C  START     COMPUTATION  BLOCK  NUMBER  31 

SM(1)=S(1) 

5M(2)=S(2) 

TTMMrTTMMPi 
C  END         COMPUTATION  BLOCK  NUMBER  31 

C  SHIFT  TO  NEW   MOVirJG  BOUNH-VRY   AND  CALC     NEW  PA^AMETFRS 

C  START     COMPUTATION  BLOCK  NUMBER  32 

IF    (      TIMEND    .GT.    TTMMPl    )    GO  TO  78 
C  END         COMPUTATION  BLOCK  NUMBER  32 

C  ==r====r==r=r=====:r=r===rr==============r===r=rr======rz=r====== 

C  START     COViPUTATION  BLOCK  NUMBER  33 

IF(JSTOP    .EQ.    1    )    GO  TO  900 
C  END         COMPUTATION  BLOCK  NUMBER  33 

MTIME=MTIME+1 
TT(MTIME)=TTMMP1 
Z0(MTIME)=70P1 

C  CALCULATE  CQNeRVED  QUANTITIES  AND  GRAPH  Op  STANDARD  S  DIMENSION 

IF(KW   .LT.      1)   GO  TO  800 
AF(1)=Y0P1 
DM=S(2) 

Af(2)=AF(1)+DM  ♦DY 
BF(1)=CP(2) 
BF(2)=C(2) 
DO  8O3  J=3»ND 
AF(J)=AF(2)+( J-2. )*DY*1G. 
K=NR  +(J-2)*10 
PF( J)=C2 (K ) 
803  CONTINUE 

DO  250  J=lrND 
XRL( J)=AF( J) 
YPL( J)=HF( J) 
250  CONTirjUE 

WRITE(6»26) 

CALL  PLOT  (ND»XPL»YRL) 
C  CONSERVATION  OF  NUMBER  CALCULATION 

Si  =C1(NL)*DZ*.5+(C1(NL)+CR(1) )*  S(l)*.5*nz 
IF(NL   .EQ.   1)    SI     =(C1(1)+CP(1) )*S(1)*.5  *DZ 
IF(NL  .GE.   2)S1  =S1  +.5*DZ*C1(1) 
IF(NL   .LE.   2)   GO  TO  221 
DO  113  J=2»NLM1 

51  =S1  +DZ*C1(J) 

113  COUTirJUE 

221     S2  =C2(riR)    *DY*.5+(C2(NR)+CB(2)  )*  S(?)*.5*DY 
DO   11*+  J=MRP1»N2M1 

52  =52  +C2(J)*DY 

114  CONTINUE 
5TM=  S1M+S2M 
ST=S1+S2 
E1=ST-STM 

RY={DF2*DT) / (DY+*2) 
E2=-DT*(C0rjPR0)  *  (CI  (l)-CO^lOUT) 

E3=RY*{C2(N2)-C2(N2Ml) )*DY 
ERR0R=E1-(F2+E3) 
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C  START     COMPUTATION  BLOCK  NUMBER  31 

800  CONTINUE 

WRlTE(6r26) 

WRITE(6f8)  TTMMPl 

WRITE(6»5)    (C1(J)    »J=1»NL»  5 

WRITE(6»41) 

WRITE(6fl)  MMfNLrNR 

WRITECeft+S) 

WRITE (6» 5)   CI (1) »CB(1) »CB(2) »DTSrDZ»DY 
WRITE(6»9)  TTMMPl 
WRITE(6»5)    (C2(J)»J=NR  »N2»20) 

NM0D=8 

NM0D1=U 

UN1=1./Z0P1 

UN2=2.*SQRT(DF2*TTMMP1) 
UN2=1./UN2 
EZ=DZ*UN1 
NLP1=NL+1 
NLP2=NLP1+1 
NLP3=NLP2+1 
AF{1)=0. 
BF(1)=C1(1) 
NC=INT(  (NL-D/NMOO) 
FNC=FLOAT(NC) 
IF(NC   .LT.     1   )   GO  TO  231 
DO  230  J=1»NC 
K=J*NMOD  +1 
JP1=J+1 
BF( JP1)=C1(K) 
AF( JP1)=(K-1)*EZ 
230  CONTINUE 

231  NCPl=NC+2 
NCP2=NC+3 
NCP3=NC+i+ 
AF(NCP1)=Z0P1*UN1 
BF(NCP1)=CB(1) 
BF(NCP2)=Y1S 
AF(NCP2)=X1S 
AF(NCP3)=0, 
BF{NCP3)=0. 

DO  252  J=1»NCP3 
XRL( J)=AF( J) 
YRL(J)=BF(J) 
252  CONTINUE 

WRlTE(6r26) 

CALL  PLOT  (NCP3rXRL»YRL) 

WRITE (6» 36) NORUN 

EY=DY*UN2 

AF(NLP2)=0. 

BF(NLP2)=CB(2) 

BF(  NLP3)=C2(NR) 

AF(NLP3)=AF(NLP2)+S(2)*EY 

DO  232  J=3»ND 

K=NLP3+J-2 

KS=NR+( J-2)*NM0ni 

BF(K)=C2(KS) 

AF (K ) =AF (NLP3) + ( J-2) *EY*NM0D1 

232  CONTINUE 

DO  234  J=1»ND 
K=J+NLP1 
AF( J)=AF(K) 
BF(J)=BF(K) 
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23**  CONTINUE 

X2SCT=X2S/2. 
NDE=0 

DO  235  J=lrND 

IF(AF(J)    .LT.     X2SCT)  NDE=NDE+1 
235  CONTINUE 

NDP1=NDE+1 
AF(NDP1)=X2S 
BF(NDP1)=Y2S 
NDP2=NDP1+1 
AF(NDP2)=  0. 
BF{NDP2)=0. 
DO  251  J=1»NDP2 
XRL( J)=AF( J) 
YRL(J)=BF(J) 
251  CONTINUE 

WRITE(6»26) 

CALL  PLOT  (N0P2»XRL»YRL) 
C  END         COMPUTATION  BLOCK  NUMBER 

C  START     COMPUTATION  BLOCK  NUMBER  35 

WRITE   (6f35)  NORUN 
PARAM{8)=CP(1) 
PARAM(9)=CB(2) 
PARAM(10)=DT*BTM 
PARAM (11) =TTMMP1*BTM 
PARAM( 12)=BDRY*BLTH 
Z0CR=Z0P1*BLTH 
WRITE(6»32) 
WRITE(6»38)  NORUN 
WRITE(6»33) 

WRITE(6»5)    (PARAM( J) » J=1.7  ) 
WRITE  (6»34) 

WRITE(6»5)  SEG» ALPHA»WFRC»  TFINaL»Z0CR  »scale? 
WRITE(6»37) 

WRITE(6»5)    {PARAM(J)    »   J=fl»12  ) 

PARAM{12)=CMAX 

PARAM( 13)=CBULK 

PARAM(14)=Y1S 

PARAM(15)=Y2S 

PARAM(16)=X25 

PARAM(17)=X1S 

PARAM(18)=BTM 

WRlTE(6r4'+) 

WRITE(6»5)    (PARAM(J)»   J=12»18  ) 
WRlTE(6»'+2) 

WRITE(6»1)   NZO»MME   »ND   »NYO»NH  tKPRINT 

C  END         COMPUTATION  BLOCK  NUMBER  35 

C  EtaRT     COMPUTATIOn'bLOCK  NUMBER  36 

J5AVE  =JSAVE  +1 

IF(  JSTOP  .EG).  JSaVE  )   GO  TO  79 
C  END         COMPUTATION  BLOCK  NUMBER  36 

J  =  JSAVE 
DF2=  DFB(J)*DFC 
BT   (J)   =(   BLTH**2)/  DF2 
WRlTE(6»5)nFB{l) »DFC»DF2»3LTH»BT( D 
TNEW(J)=  TTMMPl*  BTM/BT  (J) 
CONPRO  =  CONPR   (J)*   (   BT   (J)/  BLTH) 
CNSTA     =  CNST1(J)/BLTH 
CNSTB     =  CNST2(J)/DF2 
CNTAU     =  CNTA   (J)/BT  (J) 
DFl  =  DFA(J)/DF2 
DF2  =  l./DFC 
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WRlTE(6»5)nFB{l ) » DFC r DF? » rLTH » RT ( 1 ) »DFl,nFA(l) 
CNA  =     .5  *  CNSTA 
CNB  =  CMSTB 

CN5TX  =     SORT   (   CNA**2  +CmP*CnTAU  )-CNA 
TLIN  =  CNTAU     +  TNEW(J) 
:  CALCULATE  DT  SLOPE  ZCRIT 

DM=DZ*WFRC 

IF(   CNTAU     .LE.      .00000001)DT=   (DM*(Dm+  CnSTA)  )/CNSTR 
IF(   CNTAU     .LE.      . 0000000 1 ) DT5  =DT 
IF(   CNTAU     .LE.      .00000001)   ZCRIT  =  ZOPl 
IF(   CNTAU     .LE.      .00000001)   SLOPE  =1. 
IF(   CNTAU     .LE.      .00000001)   GO  TO  316 
ZCRIT  =  ZOPl  +  CN5TX 
SLOPE  =  CNSTX/CNTAU 
DT  =  DM/SLOPE 
DTS=DT 
316  CONTINUE 

ZONEW(J)   =  ZOPl 
TTMM=  TNEW(J) 
N0RUN=N0RUN+1 
WRITE{6»27  ) 

WRITE(6»5)TLIN»CNTAU»TTMM,   SLO^E   »ZCRIT»  CNSTV 
TW  =TNEW(J) 

TIMEND  =  TFINN(J)   *   l./BT  (J) 
WRlTE(6r39) 

WRITE(6f5)    TIMEMD»TNEW( J) » PTM » CnTAU » TLT N » BT ( J) 
ZW=  ZONEW(J) 

SCALE2  =  2.*  SQRT(DFB( J)*TFINN(J) ) 

BTM=BT( J) 

TFINAL=TFINN( J) 

CALL  TRIDNL 

PARAM(1)=DFA( J) 

PARAM(2)=DFB{ J) 

PARAM(3)=  CONPR(J) 

PARAM(U)=CONOUT 

PARAM(5)=CNST1 (J) 

PARAM(6)=CNST2( J) 

PARAM(7)=CNTA( J) 
C  ===r=====r===z=r===rrrr===r===rr=======-==r====r:r=r===: 

C  START  COMPUTATION  BLOCK  NUMBER  37 

78  CONTINUE 

IF(TIMEND   .GT.   TTMMPl)  WRITE(6»43) 

79  CONTINUE 
STOP 
END 
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11.2.     Subroutine  CENTER 


C  START     COMPUTATION  BLOCK  NUMBER  0 

IMPLICIT  DOUBLE  PRECISION  (A-H»0-Z) 
COMMON/BL20/DT»DYf DZrTHETA 

C0MM0N/BL21/DF1 »0F2 
C0MM0N/BL25/P(2) » Q ( 2 ) » PA ( 2 ) » QA ( 2 ) » C ( 2 ) » CB ( ? ) » R ( 2 ) f S ( 2 ) » SM ( ? ) 
C0MM0N/BL26/P1 (2000) »P2(2000) »Q1 (2000) » Q2 ( 2000 ) r DK 2000 )» D? ( 2000 ) 
C0MM0N/BL28/N2M1»NLM1»N2   f  NRPl»NL»MR 
COMMON/BL29/C1(2000) »C2(2000) 

COMMON  /BL30  /  CO » CK » TZERO » JEE » ALPHA » CAPA » SE6 » CMVY 
C0MM0N/RL32/C0NPR0»C0N0UT 
C  END         COMPUTATION  BLOCK  NUMBER  0 

5         FORMAT  (7Flt|.5) 

30       FORMAT  (   '+3H  S(l)   SU  SUM  FL  FR     ERR  ) 
C  CONSERVE  CALCULATION 

SUM=  DZ*   (CI (NL)+CB(1) )*.5  *5M(1)  +DY* 

1SM(2) *(C2(NR)   +CB(2) ) * . 5+C2 ( NR ) ♦DY* . 5 
C  rr=====r===============r======r=r======r===r===============r===== 

C  START     COMPUTATION  BLOCK  NUMBER  1 

P(2)=P2(NPP1) 
0(2)=Q2(NRP1) 
C  END         COMPUTATION  BLOCK  NUMBER  1 

C  ==rr=========:=r=====r=i:=============rr==rr===rr======r===rr===r= 

C  START     COMPUTATION  BLOCK  NUMBER  2 

IF(NL   .EQ.   1)   GO  TO  84° 
P(1)=P1 (NLMl) 
Q{1)=Q1(NLM1) 
C  END         COMPUTATION  BLOCK  NUMBER  2 

C  ===r======rr===rr=====r=r=rr======r=====r========r====~========= 

C  START     COMPUTATION  BLOCK   NUMBER  3 

8U9  R(l)=  2.*DT/(DZ**2) 
R (2)=2.*DF?*DT/ (DY**2) 
DO  850  L=l»2 

IF(NL   .EG.   1    .AND.   L   .EQ.   1    ) GO  TO  850 

C  COMPUTE  P  AND  Q  FOR  HALF  CELL 

DM=1.  +5(L) 

IF(L   .EO.   2   )    GO  TO  847 

DM1=DF1+S(L)/R(L) 
DM2=DF1/DM 

DM3=DM2*S(L) 

IF(L   .EO.    l)DMt|=(S(L)/R(L)  )*C1(NL) 
IF(L   .EG.   1)   GO  TO  848 
8i+7  DM1  =  1.+S(L)/R(L) 
DM2=1 ./DM 
DM3=S(L)/DM 

IF(L   .EQ.   2)DM4=(S(L)/R(L) )*C2(NR) 
848     PA(L)=DM2/(DM1-DM3*P(L) ) 

GA(L)=(DM4+Q(L)*DM3)/(DM1-DM3*P(L) ) 

850  CONTINUE 

IF(NL   .EQ.  1)PA(1)=P1(1) 
IF(NL   .EQ.  1)QA(1)=Q1(1) 
C  END         COMPUTATION  BLOCK  NUMBER  3 

C  START     COMPUTATIOn'bLOCK  NUMBER  4 

C  CALCULATION     OF  BOUNDARY  DENSITIES 

IF(NL   .EQ.   1)SAVE=  S(l) 

IF(NL   .EQ.  1)S(1)=1. 
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DM=1.+S(1) 
D11=S(1)**2 
D12=S{1)*S(2) 
RT={DZ*DF2)/(DY) 

DMl=Dll/R{  1)+  (SEG*RT*D12)/R(2)  +^^f^/n^ 

IF(NL  .EG.   1)DM1=SAVE/R(1)  +  (SEG*RT*S(2)  )/R(2.) 

DM2=-Dll/R(l)-»-(l.-S(l)  )*DF1 

IF(NL   .EQ.  1)DM2=-(SAVE/R(1)+DZ*C0NPR0) 

DM3=-RT*D12/R(2)-RT*S(1)*(1,+1,/R(2) ) 

DM4=D11*DF1/DM 

IF(NL   .EQ.  1)DM4=0. 

DM5=RT*S(1) 

DM6=S(l)*SM(l)*(Cl (NL)+C3(1) ) /R ( 1 ) +RT*SM ( 2 ) ♦S ( 1 ) *cB ( 2 ) /R ( 2 ) 
1  +S(1)*(RT*SM(2)/R(2)+RT/R(2)  )*C2(NR) 
IF(NL   .EQ.  1)DM6=DN'6+DZ*C0NPR0*C0N0UT 
IF(NL  .EQ.  1)S(1)=SAVE 
DM8=DM2+DMt+*P(l) 
DM9=DM3+DM5*P(2) 
DM10=DM6+G{1 )  *DMt++DM5*Q(2) 
DM11=DM10+QA(1)*DM8  +QA(2)*DM9 
DM12=DM1-DM8*PA ( 1 ) -SEG*DMq*P A ( 2 ) 
CB(1)=DM11/DM12 
CB(2)=SEG*CB(1 ) 
C       ^  END         COMPUTATION  BLOCK  NUMBER  U 

C  START     COMPUTATION  BLOCK  NUMBER  5 

DO  855  L=i»2 

C(L)=  PA(L)*CB(L)+QA(L) 
855  CONTINUE 

Cl(NL)=C(l) 

C2(NR)=C(2) 
C  END         COMPUTATION  BLOCK  NUMBER  5 

C  CONSERVE  CALCULATION 

SU=  DZ*S(1)*(C1(NL)+CR(1) )*.5  +   .5*  DY* 

1S(2)*(C2(NR)+CB(2) ) +C2 ( NR ) *DY* . 5 

IF{NL   .GT.l)     FL=-R(1)*(-(S{1)**2)*  Cl ( NLMl ) / ( 5 ( D +1 . ) +Cl ( ML) 
1*{S(1)-1. )+CB(l)/(S(l)+l.    )    )*DZ  *.5*DF1 

IF(NL   .EQ.1)FL=-DT  *CONPRO* (C 1 ( 1 ) -CONOUT ) 

IF(NL  .EQ.   1)FR=R(2)*   ( C2 (NRPl ) -C2 (NR ) ) *DY* . 5 

IF(NL  .EQ.  1)ERR=SU-SUM-{FL+FR) 

IF(NL   .EQ.   1)60  TO  860 

FR=R(2)*S(1)*{C2(NRP1)-C2(NR) )    *DY  *.5 
ERR=S ( 1 ) * ( SU-SUM ) - ( FL+FR ) 
860  CONTINUE 
RETURN 
END 
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.3.     Subroutine  TRIDNL 


OFORrlS  TRIDNL 

SUBROUTINE  TRIDNL 

IMPLICIT  DOUBLE  PRECISION  (A-H»0-Z) 
COMMON/RL20/DT »  DY  »  DZ »  THETa 
C0MM0N/BL21/DF1 »DF2 

C0MM0N/BL23/  B2M » B2R » B2L »  A2L » A2M » A2R »  AlL » A 1 M f  A  IP f B1 L ' RlM r R1 R 

C0MM0N/BL25/P(2) »Q(2)»PA(2)»QA(2)»C(2)»CB(2)»o(2)»S(2)»5M(2) 

C0MM0N/BL26/P1 (2000) »P2(2000) »Q1(2000) » q2 ( ?00n )» Dl ( 2000 )» D? ( 2000 ) 

C0MM0N/PL28/N2M1 f NLMl »N2    »  NRPl»NLfMR 

COMMON/BL29/C1(20  00) »C2(2000) 

C0MM0N/RL32/C0NPR0»C0N0UT 
5         FORMAT  (TFll.S) 
C  CALC.   P»S  AND  DUMMIES 

C  START     COMPUTATION  BLOCK  NUMBER  1 

CALL  INTER 

C  END         COMPUTATION  BLOCK  NUMBER  1 

IF(NL   .EQ.  1)0M1=S(1)*DZ 

IF(NL   .EG.  1)OM2=Dm1*DM1/(2,*DT) 

IF(NL   .EQ.   l)DM=i./(  DF1+DM2+0M1*C0MPR0) 

IFCNL   .EQ.   1    )P1  (1)=DM*D'='1 
C  START     COMPUTATION  BLOCK  NUMBER  2 

IF(NL   .GT.   1   )DM=1./(DZ*DZ/(2.*DT)   +D7*  CONPRO+DFl  ) 

IFCNL   .GT.   1    )P1  (DrDFl+n-^ 
C  r============r=rr=====r======r=r========r=r==r=========rr=r======= 

C  END         COMPUTATION  BLOCK  NUMBER  2 

C  START     COMPUTATIOrj  BLOCK  NUMBER  3 

IF(NL   .LE.   2   )G0  TO  106 

DO  lOS  J=  2fNLMl 

JM1=J-1 

DKJ)   =A2M+P1  ( JM1)*A2L 
PI ( J)=-A2R/D1 ( J) 

105  CONTIrjUE 

106  CONTINUE 

C  END         COMPUTATION  BLOCK  NUMBER  3 

C  START     COMPUTATION  BLOCK  NUMBER  4 

DO  101   J=  1»N2M1 
K=rJ2-J 
KP1=K+1 

D2(K)=R2M+P2(KP1) 
P2(K)=-1./D2(K) 
101  CONTINUE 
C  END         COMPUTATION  BLOCK  NUMBER  4 

RETURN 
END 
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Subroutine  INTER 


SUBROUTINE  INTER 

IMPLICIT  DOUBLE  PRECISION  (A-H»0-Z) 
IMPLICIT  DOUBLE  PRECISION  (A-H»0-Z) 
COMMON/BL20/DT»DY»DZ»THETA 
C0MM0N/BL21/DF1»DF2 

C0MM0N/BL23/  B2M r B2R r B2L f A2L » A2M » A2R r A IL » AIM » AlR » BIL » RIM » BIR 
5         FORMAT  (TFU.S) 
RA1=DT/(DZ**2) 
RA2=DT/(DY**2) 
A2L=  DFl 
A2R=A2L 

A2M=-(2.*DF1+1./RA1) 

A1L=1. 

A1R=  AIL 

A1M=-1 ,/RAl 

B2L=1. 

B2R=B2L 

B2M=-(2.+l./(DF2*RA2)  ) 
B1L=1. 

BIR  =B1L 

BIM  =  -l./(DF2*RA2) 

RETURN 
END 
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12 .  GLOSSARY 


See  eqs   (5.6)   and  (5.7)   and  table  5. 

An  array  that  stores  the  inhomogeneous  terms  in  the  tridiago- 
nal  equations  for  the  z-region   (table  5) . 

Swelling  ratio  a   [eq  (2.5)]. 
See  eqs   (5.10)   and  table  7. 

The  point  representing  infinity  in  the  y-region   [eq  (2.11^)]. 

The  unit  being  used  to  dimensionalize  lengths    (sec.  4.3). 

The  unit  being  used  to  dimensionalize  time   (sec.  4.3). 

An  array  that  stores  the  inhomogeneous  terms  in  tridiagonal 
equations  in  the  y-region   [eq  (5.10)   and  table  7]. 

Material  parameters   [eqs   (2.11),    (2.12),  and  (2.13)]. 

Arrays  representing  concentrations  at  the  grid  points  in  the 
z-  and  y-regions,  respectively.     At  the  beginning  of  an  itera- 
tion of  loop  78,  CI (  )   represents  Ci (   ,  TTMM),  and  at  the  end 
of  an  iteration,  CI (  )   represents  C^ (   ,  TTMMPl) .  Likewise, 
C2(  )  . 

These  represent  C^ (zq (TTMMPl) ,  TTMMPl),  i  =  1,  2,  i.e.,  the 
concentrations  in  the  oxide  and  silicon  at  the  moving  boun- 
dary. 

The  concentration  in  the  bulk  of  the  silicon  and  at  y  =  0 
[eqs   (2.13)   and   (2.19),  Cg  =  CBULK] . 

A  subroutine  for  calculating  CI (NL) ,   C2 (NR)  ,  CB(1),  CB(2) 
(sees.   7.2  and  8.2). 

The  value  of  C2(y,0)   for  y  £  NH  in  eq  (2.13). 
By  definition,  A/2. 
Another  name  for  CNSTB. 

The  constants  A  and  B,  respectively,  in  eqs   (4.5)   to  (4.9) 
(table  2) . 

Defined  in  eq   (4.9)   and  FORTRAN  name  for  ZCRIT. 
The  constant  t  in  eqs   (4.8)   and  (4.9). 
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CONOUT 


CONPRO 

Dl 

D2 

DFl 
DF2 
DFC 

DM1  I 
through  > 
DM4  ' 

DM1  1 

through  > 
DM10  J 

DT 

DTS 

DY 

DZ 

EY  I 
EZ  ' 

FNL 

FNR 

FNZO 

JSAVE 

JSTOP 
KPOW 

KPRINT 

LOOP  78 
m 


-  A  constant  in  eq  (2.18)   representing  the  boron  concentration 
in  the  oxygen  at  the  oxygen-oxide  interface   (sec.  5.7); 
CONOUT  =  Cout- 

-  A  proportionality  constant  in  eq  (2.18)  representing  the 
evaporation  rate  of  boron  from  the  oxide  into  the  oxygen 
ambient  (sec.  5.7);  CONPRO  =  Cp. 

-  An  intermediate  array  used  to  compute  the  arrays  Ql ,  PI 
[eq   (6.9')   and  sec.   7.3,   computational  block  #3]. 

-  An  intermediate  array  used  to  compute  arrays  P2,  Q2  [sec. 
7.3,  block  #4,  and  eq  (6.16')]. 

-  The  diffusion  coefficient  Di  in  the  oxide  or  z-region. 

-  The  diffusion  coefficient  D2  in  the  silicon  or  y-region. 

-  A  normalizing  factor  set  equal  to  one  in  VERSION  1. 

-  For  these  quantities  used  in  computational  block  #3  of  siib- 
routine  CENTER,   see  eqs   (6.I0M   and  (6.17'). 

-  For  these  quantities  used  in  computational  block  #4  of  sub- 
routine CENTER,   see  eqs    (6.20),    (6.21),   and  (6.22). 

-  The  time  mesh  width  At    (sec.   4.4) . 

-  A  value  of  DT  that  is  stored  before  the  moving  boundary  zq 
crosses  a  grid  point   (sec.   7.1,  computational  block  #20). 

-  The  mesh  width  Ay  in  the  y-region   (sec.  4.1) . 

-  The  mesh  width  Az  in  the  z-region    (sec.  4.1). 

-  Scaled  value  of  DY,   DZ  used  in  plot  output. 

-  The  real  value  of  NL  -  1. 

-  The  real  value  of  NR  -  1. 

-  The  quantity  NZO  floated,  i.e.,  made  real. 

-  A  parameter  used  to  shift  control  in  the  program  (computa- 
tional block  #36,    sec.    7.1).     JSAVE  =  0  in  VERSION  1  to 
begin  with. 

-  In  VERSION  1,   JSTOP  equals  1. 

-  A  signal  to  increase  the  value  of  DT  when  KPOW  =  100  and 
ZOACEL  <_  0.2    (computational  block  #16,   sec.  7.1). 

-  The  concentrations  and  certain  other  integer  parameters  are 
printed  out  when  MM  =  0  modulo    (KPRINT) . 

-  Each  iteration  of  this  loop  calculates  the  concentrations  at 
one  time.   At.      (Later  see  computational  block  #14,   sec.  7.1). 

-  Another  symbol  for  segregation  coefficient  SEG   [eq  (2.2)]. 
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This  is  the  "DO"  index  in  loop  78. 

This  integer  specifies  the  number  of  iterations  in  loop  78 
(computational  block  #14,   sec.  7.1). 

The  value  of  NL  at  time  t  =  0  in  VERSION  1  (computational 
block  #9,   sec.   7.1) . 

The  number  of  fixed  grid  points  in  the  y-region,  N2  =  NYO  + 
1   [computational  block  #9,  sec.  7.1,  and  eq  (4.2)]. 

Tells  how  many  y-grid  points  to  the  right  of  NR  will  be 
plotted  in  output  graph   (block  #34,   sec.  7.1). 

A  loop  limit  index  in  block  #26,   sec.  7.1. 

The  point  on  the  y-axis  where  the  initial  concentrations  in 
the  silicon  changes  from  CMAX  to  CBULK   [eq   (2.13)   and  sec. 
7.1,  computational  block  #10]. 

See  end  of  section  4.1. 

The  z-grid  points  are  plotted  modulo  NMOD  (block  #34,  sec. 
7.1)  . 

The  y-grid  points  from  NR  to  the  right  are  plotted  modulo 
NMODl   (block  #34,   sec.  7.1). 

In  VERSION  1.     NORUN  =  1. 

See  end  of  section  4.1. 

The  number  of  mesh  widths  of  length  DY  that  cover   (0,  BDRY) 
(sec.  4.1)  . 

A  parameter  used  in  the  definition  of  DZ    [eq  (4.3)]. 

These  are  the  arrays  whose  functions  are  described  in  eqs 
(3.5)   and   (3.7).     The  quantities  with  the  exception  of  Pl(l) 
are  computed  in  subroutine  TRIDNL.     See  eqs   (6.5'),  (6.5'), 
and  (6.7')   for  Pl(l)   and   (6.9')   for  P1(J),  J  >  1.     See  com- 
putational blocks  #23  and  #24  for  computation  of  Pl(l). 
See  eq  (6.16')   for  the  derivation  of  (6.16). 

Alternate  names  for  Pl(NLMl),  P2  (NRPl)    [see  eq  (6.18')]. 

Coefficient  linking  CI (NL)   with  CB(1)   and  C2 (NR)   with  CB(2) 
calculated  in  subroutine  CENTER,  block  #3    [eqs   (6.19')  and 

(5.12')].     Function  of  this  array  is  explained  in  section 

3.1,  steps  #2  and  #3. 

An  array  containing  various  material  input  constants  and 
computer  program  constants  of  possible  interest  that  are 
program  output   (computational  blocks  #2  and  #35,  sec.  7.1. 

An  NBS  "in-house"  plotting  program  available  from  NBS  Applied 
Mathematics  Division.     Documented  in  Appendix  2. 

Everything  to  do  with  these  arrays  is  explained  in  the  same 
places  as  Pi  and  P2 . 

Ill 


Q(l)l 
Q(2)» 

R 


R(l)) 
R(2)> 

RAl) 
RA2' 


RT 


S(l)i 
S(2)! 


SAVE 

SCALE2 

SEG 

SES 

SLOPE 

SM(1)) 
SM(2) ' 

TFINAL 
TIMEND 

TJ 

TLIN 
TRIDNL 
TTMM 
TTMMPl 


-  Alternate  names  for  Ql(NLMl),  Q2 (NRPl)    [eq  (6.18')]. 

-  Another  name  for  RT. 

-  These  are  quantities  that  are  defined  in  eqs   (5.35)  and 
calculated  in  computational  block  #3  of  subroutine  CENTER. 
[See  also  eqs   (5.25)   and   (5.16,  5.16')]. 

-  Combinations  of  mesh  widths  used  in  subroutine  INTER  (tables 
5  and  7) . 

-  This  combination  of  mesh  quantities  is  defined  in  eq  (5.35) 
and  appears  also  in  eqs   (5.38)   and  (5.34)  (computational 
block  #4  of  subroutine  CENTER) . 

-  See  eqs   (5.35),    (5.29),    (5.22),    (5.15)   and  figures  5  and  6. 
These  quantities  are  computed  in  computational  block  #21, 
section  7.1. 

-  See  computational  block  #4  of  subroutine  CENTER;   see  remarks 
at  the  end  of  section  5.9  before  the  SUMMARY. 

-  Essentially  the  same  as  UN2 ,  a  scale  factor  used  in  plotting. 
(See  sec.  4.4) . 

-  The  FORTRAN  symbol  for  m,  the  segregation  coefficient  [eq 
(2.2) ]  . 

-  A  signal  which  when  negative  implies  the  moving  boundary  yo 
reached  a  new  y-grid  point  during  the  previous  iteration  of 
loop  78   (computational  block  #19,   sec.  7.1). 

-  The  speed  of  the  moving  boundary  in  dry  oxidations    [eq  (4.8)] 

-  See  the  same  equations  and  figures  cited  in  description  of 
S(l),  S (2) .     These  quantities  are  computed  in  computational 
block  #12  and  #31,   section  7.1. 

-  This  is  the  time  in  hours  that  is  inputted  into  the  program. 

-  When  TTMMPl  is  greater  than  TIMEND   (which  is  TFINAL  in  BTM 
units) ,  the  program  stops  and  control  is  shifted  to  printing 
output   (computational  blocks  #32  and  #4,   sec.  7.1). 

-  A  signal  which,  when  its  value  is  greater  than  or  equal  to 
1 ,  indicates  that  the  moving  boundary  Zq  has  reached  a  new 
z-grid  point  during  the  previous  iteration  of  loop  78-. 

-  The  quantity  T  in  dry  oxidations   (table  4.1);   also  called 
CNTAU. 

-  The  subroutine  where  arrays  PI,  Dl  are  calculated  (sees. 
7.3  and  8.3) . 

-  The  time  at  the  beginning  of  the  MMth  iteration  of  loop  78 
in  units  of  BTM  (computational  block  #18,    sec.  7.1). 

-  The  time  at  the  end  of  the  MMth  iteration  of  loop  78   (in  BTM 
units) . 
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This  quantity  is  always  zero  in  VERSION  1. 

Scaling  factors  used  in  the  plot  output   (computational  block 
#34,   sec.    7.1,   and  sec.  4.4). 

The  fraction  of  DZ  that  the  moving  boundary  Zq  moves  in  the 
first  iteration  of  loop  78   [eq  (4.4)   and  computational  block 
#7,   sec.   7.1]  . 

Quantities  used  in  setting  scales  for  plots  of  concentrations 
(computational  block  #34,   sec.   7.1,  and  sec.  4.4). 

A  cutoff  for  points  with  large  abscissas   (sec.  4.4). 

An  array  for  abscissas  used  in  plotting   (sec.  3.4). 

The  position  at  time  t  in  the  y-coordinate  frame,  calculated 
by  using  eq  (2.8). 

Quantities  used  in  setting  scales  for  plots  of  concentra- 
tions  (computational  block  #34,   sec.   7,1,  and  sec.  4.4). 

Always  zero  in  VERSION  1. 

Position  in  y-region  of  moving  boundary  at  time,  TTMMPl  (com- 
putational block  #18,   sec.  7.1). 

An  array  for  ordinates  used  in  plotting   (sec.  3.4). 

The  position  of  the  moving  boundary  in  the  z-coordinate 
frame  at  time  t   (sec.  4.2). 

Same  as  CNSTX. 

Always  zero  in  VERSION  1. 

When  the  moving  boundary  is  moving  slowly,  this  parameter  is 
small  and  the  value  of  DT  is  increased  in  order  to  speed  up 
distance  the  boundary  moves  in  a  given  iteration  of  loop  78 
(computational  block  #16,   sec.  7.1). 

ZOPl  in  units  of  micrometers   (computational  block  #34,  sec. 
7.1)  . 

Position  in  z-region  of  moving  boundary  at  time,  TTMMPl  (com- 
putational block  #18,   sec.  7.1). 

See  computational  block  #20,  sec.  7.1. 

See  computational  block  #16,  sec.  7.1. 

Always  zero  in  VERSION  1. 
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Appendix  1 


FORMULAS  USED  IN  DERIVATION  OF  DISCRETE  ALGEBRAIC  EQUATIONS 

We  collect  here  for  convenience  some  standard  finite  difference  approxi- 
mations used  in  the  discretizations  effected  in  chapter  five.     The  deriva 
tions  of  these  formulas  can  be  found  in  a  text  book  on  numerical  analysis 


1.     TAYLOR'S  SERIES 

f  (x+Ax)   =  f  (x)   +  f '  (x)  Ax  +   Ax2  +.  .  .+   Ax"  +.  .  .  (Al) 

2  nl 


2.     QUADRATURE  APPROXIMATIONS 


a.     Midpoint  rule 


f  (x)dx  = 


f(^) 


(b-a)  +  0[(b-a)3] 


(A2) 


b.     Right  end  point 


f(x)dx  =  f(b)Ax  +  0[(b-a)2] 


(A3) 


c.     Left  end  point 

/  f(x)dx  =  f(a)  (b-a)  +  0[(b-a)2]  (A4) 
a 


d.     Trapezoid  rule 


f(x)dx  =  ^^^^^^^^^    (b-a)   +  0[(b-a)3]  (a5) 


3.     DERIVATIVE  APPROXIMATIONS 


a.     Forward  differences 


df (x)   ^  f(x+Ax)   -  f(x) 
dx  Ax 


+  0[Ax] 


(A6) 
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Backward  differences 
df(x)       f(x)   -  f(x-Ax) 


dx  Ax 

Central  differences 


+  0[Ax]  (A?) 


df  (x)  f  (x+Ax)-f  (x-Ax)  ^  2i 
— :;          -   —  +  0[Ax^] 


-,  _  .  -   .  (A8) 

dx  2Ax 

4.     THREE-POINT  APPROXIMATIONS 

Let  X2 ,  X,  X]^  be  three  points  such  that  xj  >  x  >  X2 ;  then  a  three-point 
approximation  to  the  derivative  of  t(x)   at  point  x  is 

,^          Ax.   f(x-Ax„)          (Ax  -Ax^)f(x) 
df  _  _       1  2      _^         1       2   ^ 

dx  Ax^   (Ax^+Ax^)  Ax^  Ax^ 

(A9) 

Ax^  f (x-Ax^) 


Ax^  (Ax^+Ax^) 


where  Ax^  =  x^^  -  x  and  Ax^  =  x  -  x^ 
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APPENDIX  2 


Plotting  vuh routines 
B.  L.  Joiner  and  S.  T.  Peavy 


Five  FORTRAN  subroutines  for  producing  plots  similar  to  the 
one  shoun  in  Figjre  1  are  now  available  on  the  NBS  UNI VAC  1108 
FASTRAND  file  named  PLOTS. 

Questions  related  to  the  operation  of  these  subroutines  may- 
be directed  to: 

Brian  L.  Joiner  or  Sally  T.  Peavy 
A337  Administration  Bldg. 
National  Bureau  of  Standards 
Washington,  D.  C.  20234 
921-2315 

Task  numbers  for  the  IINIVAC  1108  may  be  obtained  from  the 

Ccmputer  Services  Division 
A238  Adninistration  Bldg. 
National  Bureau  of  Standards 
Washington,  D.  C.  20234 
921-3364 

General  Remarks 

Except  as  specified  below,  the  suLx-outines  automatically 
figure  out  limits  for  the  plots  based  on  the  smallest  and  largest 
data  points.    A  new  page  is  not  called  by  any  of  these 
subroutines.    This  allows  the  user  to  label  the  top  of  the  page 
before  calling  for  tl^e  plot.    These  subroutines  do  not  change 
any  of  the  values  of  tlie  arguments.    Plotting  is  done  by 
repetitively  seardung  the  arrays  rather  than  by  sorting  them. 

The  resolution  of  the  plots  is  51  characters  high  by  101 
characters  v.ide  and  each  plot  consumes  54  lines  counting  borders 
and  scale  labeling.    The  length  and  v,ldth  of  the  actual  plotting 
area  are  Sh  inches  and  10  inches  respectively,  and  the  overall 
dimensions  including  borders  and  scale  labeling  are  9  inches 
and  11'4  inches  respectively. 
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These  five  subroutines  may  be  called  in  from  FASTRAND  by  inserting 
the  follov/ing  pair  of  cards  after  the  RUN  card: 

e    XQT  CUR 
INF  PLOTS 

viiere  @  means  a  7  and  8  punched  in  card  column  one. 


Subroutines 

I.    PLOT  (N,  X,  Y) 

Plots  the  data  in  one  column  versus  that  in  another  column. 

N       :  -The  number  of  points  to  be  plotted. 

X       :    A  vector  (one  dimensional  array)  containing  the 
values  of  the  abscissa. 

Y       :    A  vector  containing  the  values  of  the  ordinate. 

Example  of  main  program. 

DIMENSION    X(10),  Y(10) 
DO    2    I  =  1,5 
X(I)  =  I 
2      Y(I)  =  SIN(X(I)) 
N  =  5 

CALL  PLOT  CN,  X,  Y) 
STOP 


II.    PLOTS  G^^GS,  X,  Y,  NRMX,  NROW) 

Plots  the  data  in  up  to  5  pairs  of  columns.    The  symbols 
used  are    .  *  +  ,  -    respectively  for  the  5  curves. 

NARGS    :    The  number  of  curves  to  be  plotted. 

X  :    A  matrix  (t\co  dimensional  array)  having  up  to 

5  coluTjis  and  containing  the  values  of  the 
abscissa. 

Y  :    A  matrix  having  the  sa-ne  dimensions  as  X  and 

containing  the  values  of  the  ordinate.  The 
first  colunm  of  Y  is  plotted  versus  the  first 
column  of  X,  etc. 

NRMX     :    A  vector  (one  dimensional  array)  containing  the 
number  of  points  to  be  plotted  in  each  of  the 
(5  or  less)  column  pairs. 

NROW     :    The  number  of  rows  specified  for  X  and  Y  in  their 
dimension  statement (s)  in  the  main  program. 

(An  example  is  given  on  the  following  page.) 
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Example  of  main  program. 


DINENSION   XC85,4),  YC85,4),  NRNKC4) 

DO    5    I  =  1,20 

X(I,1)  =  I-IO 
-5      Y(I,1)  =  ABS(X(I,1)) 

DO    10    I  =  1,37 

X(I,2)  =  1-20 
10     YCI,2)  =  20.  -  X(I,2)**2 

NARGS  =  2 

NR^IX  (1)  =  20 

NRMX  (2)  =  37 

NROW  =85 

CALL  PLOTS  OIARGS,  X,  Y,  NRMX,  NRO'.V) 
STOP 

III.    PLOTL  (N,  X,  Y,  IT) 

Plots  the  data  in  one  column  versus  that  in  another  column 
using  the  symbols  specified  in  the  IT  column. 

N       :    The  number  of  points  to  be  plotted. 

X       :    A  vector  (one  dimensional  array)  containing 
the  values  of  the  abscissa. 

Y       :    A  vector  containing  the  values  of  the  ordinate. 

IT     :    A  vector  containing  numbers  bet\veen  1  and  26 
which  dictate  the  letter  of  the  alphabet  to 
be  used  as  a  plotting  symbol. 

number:  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 
letter:    ABCDEFGHI     J     K     L     M     N  0 

16    17    18    19    20    21    22    23    24    25  26 
PQRSTUVWXYZ 

Example  of  main  program. 

DIMENSION    X(IOO),  Y(IOO),  ITClOO) 
DO    20    I  =  1,5 
X(I)  =  I 
Y(I)  =  SIN(XCI)) 
20     IT(I)  =  19 

DO    30    I  =  6,10 
X(I)  =  1-5 
Y(I)  =  coscxci)) 

30      ITCD  =  3 

CALL  PLOTL  (10,  X,  Y,  IT) 
STOP 
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IV.    PLOTIA  CN,  X,  Y,  IT,  YMINS,  YMAXS) 


This  subroutine  is  the  same  as  PLOTL  except  that  the 
limits  for  the  Y  axis  (ordinate)  are  specified  in  the  call 
Statement.    If  any  points  fall  outside  the  specified  limits, 
the  limits  are  stretched  so  as  to  include  all  points. 

Example  of  main  program  is  same  as  for  PLOTL  except  for  call. 

CALL  PLOTLA  (10,  X,  Y,  IT,  -1.,  1.) 


V.    PLOTLF  (N,  X,  Y,  IT,  XMIN,  XMOX,  YMIN,  YMAX) 

This  subroutine  is  also  similar  to  PLOTL  except  that  here 
the  limits  for  both  X  and  Y  are  specified  in  the  call  statement. 
Any  points  falling  outside. the  specified  limits  are  omitted  but 
a  talley  is  kept  and  the  number  of  offending  points  is  printed 
below  the  plot. 

Exaiiple  of  main  program  is  same  as  for  PLOTL  except  for  call. 
CALL  PLOTLF  (10,  X,  Y,  IT,  0.,  5.,  -1.,  1.) 
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|:omputer  sciences.  Papers  cover  a  broad  range  of  subjects, 
(With  major  emphasis  on  measurement  methodology,  and 
lie  basic  technology  underlying  standardization.  Also  in- 
j:Iuded  from  time  to  time  are  survey  articles  on  topics  closely 
(-elated  to  the  Bureau's  technical  and  scientific  programs.  As 
^  special  service  to  subscribers  each  issue  contains  complete 
|dtations  to  all  recent  NBS  publications  in  NBS  and  non- 
NBS  media.  Issued  six  times  a  year.  Annual  subscription: 
{domestic  $17.00;  foreign  $21.25.  Single  copy,  $3.00  domestic; 
53.75  foreign. 

Note:  The  Journal  was  formerly  published  in  two  sections: 
Section  A  "Physics  and  Chemistry"  and  Section  B  "Mathe- 
matical Sciences." 
>DI  MEN  SIGNS/ NBS 

This  monthly  magazine  is  published  to  inform  scientists, 
engineers,  businessmen,  industry,  teachers,  students,  and 
jconsumers  of  the  latest  advances  in  science  and  technology, 
iwith  primary  emphasis  on  the  work  at  NBS.  The  magazine 
highlights  and  reviews  such  issues  as  energy  research,  fire 
iprotection,  building  technology,  metric  conversion,  pollution 
labatement,  health  and  safety,  and  consumer  product  per- 
iformance.  In  addition,  it  reports  the  results  of  Bureau  pro- 
grams in  measurement  standards  and  techniques,  properties 
.of  matter  and  materials,  engineering  standards  and  services, 
linstrumentation,  and  automatic  data  processing. 

Annual  subscription:   Domestic,  $1 1.00;  Foreign  $13.75 

'  NONPERIODICALS 

(Monographs — Major  contributions  to  the  technical  liter- 
ature on  various  subjects  related  to  the  Bureau's  scientific 
and  technical  activities. 

iHandbooks — Recommended  codes  of  engineering  and  indus- 
tri:il  practice  (including  safety  codes)  developed  in  coopera- 

1  with  interested  industries,  professional  organizations, 

J  regulatory  bodies. 
Special  Publications — Include  proceedings  of  conferences 
^p>)nso^ed  by  NBS,  NBS  annual  reports,  and  other  special 
publications  appropriate  to  this  grouping  such  as  wall  charts, 
pocket  cards,  and  bibliographies. 

Applied  Mathematics  Series — Mathematical  tables,  man- 
uals, and  studies  of  special  interest  to  physicists,  engineers, 
chemists,  biologists,  mathematicians,  computer  programmers, 
and  others  engaged  in  scientific  and  technical  work. 
National  Standard  Reference  Data  Series — Provides  quanti- 
tative data  on  the  physical  and  chemical  properties  of 
materials,  compiled  from  the  world's  literature  and  critically 
evaluated.  Developed  under  a  world-wide  program  co- 
ordinated by  NBS.  Program  under  authority  of  National 
Standard  Data  Act  (Public  Law  90-396). 


NOTE:  At  present  the  principal  publication  outlet  for  these 
data  is  the  Journal  of  Physical  and  Chemical  Reference 
Data  (JPCRD)  published  quarterly  for  NBS  by  the  Ameri- 
can Chemical  Society  (ACS)  and  the  American  Institute  of 
Physics  (AIP).  Subscriptions,  reprints,  and  supplements 
available  from  ACS,  1155  Sixteenth  St.  N.W.,  Wash.,  D.C. 
20056. 

Building  Science  Series — Disseminates  technical  information 
developed  at  the  Bureau  on  building  materials,  components, 
systems,  and  whole  structures.  The  series  presents  research 
results,  test  methods,  and  performance  criteria  related  to  the 
structural  and  environmental  functions  and  the  durability 
and  safety  characteristics  of  building  elements  and  systems. 
Technical  Notes — Studies  or  reports  which  are  complete  in 
themselves  but  restrictive  in  their  treatment  of  a  subject. 
Analogous  to  monographs  but  not  so  comprehensive  in 
scope  or  definitive  in  treatment  of  the  subject  area.  Often 
serve  as  a  vehicle  for  final  reports  of  work  performed  at 
NBS  under  the  sponsorship  of  other  government  agencies. 
Voluntary  Product  Standards — Developed  under  procedures 
published  by  the  Department  of  Commerce  in  Part  10, 
Title  15,  of  the  Code  of  Federal  Regulations.  The  purpose 
of  the  standards  is  to  establish  nationally  recognized  require- 
ments for  products,  and  to  provide  all  concerned  interests 
with  a  basis  for  common  understanding  of  the  characteristics 
of  the  products.  NBS  administers  this  program  as  a  supple- 
ment to  the  activities  of  the  private  sector  standardizing 
organizations. 

Consumer  Information  Series — Practical  information,  based 
on  NBS  research  and  experience,  covering  areas  of  interest 
to  the  consumer.  Easily  understandable  language  and 
illustrations  provide  useful  background  knowledge  for  shop- 
ping in  today's  technological  marketplace. 
Order  above  NBS  publications  from:  Superintendent  of 
Documents,  Government  Printing  Office,  Washington,  D.C. 
20402. 

Order  following  NBS  publications— NBSlR's  and  FIPS  from 
the  National  Technical  Information  Services,  Springfield, 
Va.  22161. 

Federal  Information  Processing  Standards  Publications 
(FIPS  PUB) — Publications  in  this  series  collectively  consti- 
tute the  Federal  Information  Processing  Standards  Register. 
Register  serves  as  the  official  source  of  information  in  the 
Federal  Government  regarding  standards  issued  by  NBS 
pursuant  to  the  Federal  Property  and  Administrative  Serv- 
ices Act  of  1949  as  amended.  Public  Law  89-306  (79  Stat. 
1127),  and  as  implemented  by  Executive  Order  11717 
(38  FR  12315,  dated  May  11,  1973)  and  Part  6  of  TiUe  15 
CFR  (Code  of  Federal  Regulations). 

NBS  Interagency  Reports  (NBSIR) — A  special  series  of 
interim  or  final  reports  on  work  performed  by  NBS  for 
outside  sponsors  (both  government  and  non-government). 
In  general,  initial  distribution  is  handled  by  the  sponsor; 
public  distribution  is  by  the  National  Technical  Information 
Services  (Springfield,  Va.  22161)  in  paper  copy  or  microfiche 
form. 


BIBLIOGRAPHIC  SUBSCRIPTION  SERVICES 


The  following  current-awareness  and  literature-survey  bibli- 
ographies are  issued  periodically  by  the  Bureau: 
Cryogenic  Data  Center  Current  Awareness  Service.  A  litera- 
ture survey  issued  biweekly.  Annual  subscription:  Domes- 
tic, $25.00;  Foreign,  $30.00. 
Liquified  Natural  Gas.  A  literature  survey  issued  quarterly. 
Annual  subscription:  $20.00. 


Superconducting  Devices  and  Materials.  A  literature  sirrvey 
issued  quarterly.  Annual  subscription:  $30.00.  Send  subscrip- 
tion orders  and  remittances  for  the  preceding  bibliographic 
services  to  National  Bureau  of  Standards,  Cryogenic  Data 
Center  (275.02)  Boulder,  Colorado  80302. 
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